Monte Carlo Simulations as a Decision Support to Interpret δ15N Values of Nitrate in Groundwater

Intense farming is often associated with the excessive use of manure or fertilizers and the subsequent deterioration of the groundwater quality in many aquifers worldwide. Stable isotopes of dissolved nitrate (δ15N and δ18O) are widely used to determine sources of nitrate contamination and denitrification processes in groundwater but are often difficult to interpret. Thus, Monte Carlo simulations were carried out for a site in lower Bavaria, Germany, in order to explain δ15N observations in a porous groundwater system with two aquifers, the main aquifer (MA) and several smaller perched aquifers (PA). For evaluating potential contributions, frequency distributions of δ15N were simulated deriving from (I) the mixing of different nitrate sources, related to land use, as input to groundwater, combined with (II) transport of nitrate in groundwater and (III) microbial denitrification. Simulation results indicate a source‐driven isotopic shift to heavier δ15N values of nitrate in groundwater, which may be explained by land use changes toward a more intensified agriculture releasing high amounts of manure. Microbial denitrification may play a role in the PA, with simulated δ15N distributions close to the observations. Denitrification processes are however unlikely for the MA, as reasonable simulation curve fits for such a scenario were obtained predominantly for unrealistic portions of nitrate sources and related land use. The applied approach can be used to qualitatively and quantitatively evaluate the influence of different potential contributions, which might mask each other due to overlapping δ15N ranges, and it can support the estimation of nitrate input related to land use.


Introduction
Many countries around the world rely on groundwater for safe drinking water supply. With the increasing pressures of climate change, population growth and rapidly developing economies on the quality and quantity of water resources, groundwater protection is a global challenge. Therefore, groundwater monitoring as well as efficient water saving strategies and pollution management must be viewed along with the analysis of the selfpurification potential.
In Europe, several regulatory frameworks aim to improve surface water and groundwater quality (European Commission 2000, 2015. Nevertheless, nitrate is still one of the main non-point source contaminants in groundwater and concentrations are widely elevated in the United States, Germany, Denmark, and elsewhere (Jørgensen and Stockmarr 2009;Hansen et al. 2012;Werner and O'Doherty 2012). In Germany, 28% of regularly sampled groundwater wells (EU monitoring network) showed nitrate concentrations exceeding the European drinking water limit of 50 mg/L in the years 2012 to 2015 (Keppner et al. 2016).
Sorption of contaminants to organic matter or mineral surfaces and microbial degradation in groundwater have the potential to significantly reduce the concentrations of pollutants in groundwater (Archna and Sobti 2012;Clark 2015). However, such processes are strongly dependent of the aquifer's properties, the microbial community and the redox conditions in the aquifer (Hiscock 2005). In a recently published study, it was shown that the degradation potential of redox sensitive pollutants and the nitrate vulnerability of aquifers may be best assessed by the determination of O 2 reduction rates and denitrification lag times in groundwater (Wild et al. 2018). However, if the availability of such data is limited NGWA.org Vol. 58, No. 4-Groundwater-July-August 2020 (pages 571-582) for a catchment area, it is difficult to characterize the selfpurification potential of the groundwater system. Under such situations, the microbial reduction of nitrate can be evaluated using stable isotope ratios of δ 15 N and δ 18 O of nitrate and may enable the identification of aquifers that are more susceptible to nitrate contamination (Aravena et al. 1993;Kendall and McDonnell 1998;Böhlke et al. 2002;Einsiedl and Mayer 2006;Wild et al. 2018). In several field studies, it was demonstrated that during denitrification both δ 15 N and δ 18 O of nitrate increase, that is, enrich in heavier isotopes, while nitrate concentrations decrease (Kendall and McDonnell 1998). However, there have also been ambiguities and difficulties interpreting δ 15 N and δ 18 O values of nitrate in groundwater systems. For instance, the mixing of different unreacted nitrate sources such as manure (elevated δ 15 N > 7‰, low δ 18 O ≤ 5‰) with unreacted nitrate deriving from precipitation (low δ 15 N ∼0‰, elevated δ 18 O ∼60‰) (Kendall and McDonnell 1998) can be misleadingly interpreted as microbial denitrification (Pauwels et al. 2000;Xue et al. 2009). Therefore, the identification of denitrification via a characteristic slope in a 2D isotope plot (δ 18 O versus δ 15 N) may often fall short for groundwater systems impacted by a mixture of different nitrate sources. In addition, during denitrification, δ 18 O originating from ambient water may be incorporated into dissolved nitrate by back reactions of NO 2 − to NO 3 − and can overprint the expected enrichment of 18 O in the remaining nitrate as reported in literature (Wunderlich et al. 2013;Granger and Wankel 2016). This implies that there is no typical slope as a robust diagnostic tool for the characterization of denitrification under environmental conditions as often suggested in literature (e.g., Amberger and Schmidt 1987;Boettcher et al. 1990). Consequently, we hypothesize that statistical analysis can help as an additional tool for the interpretation of δ 15 N of nitrate in groundwater. To reach this goal, we assessed possible influences on δ 15 N values including the mixing of different nitrate sources (each with characteristic δ 15 N values), linked with hydrodynamic processes and microbial denitrification.
Transport modeling often applies numerical solutions to include heterogeneities (Cirpka and Helmig 1999). However, in many study areas calibration may be difficult due to a low spatial resolution of known aquifer properties and details about geology. Otherwise, numerical modeling requires a stochastic framework for uncertainty analysis (Simmons et al. 1995). Literature also shows that if an extended data set for a groundwater system is missing, it makes sense to use simple lumped-parameter models Zuber 1982, 1996), which also use a statistical characterization of the variability of groundwater ages, notwithstanding the many other sources of uncertainty. Isotope mixing models implementing a Bayesian framework are widely used in ecological food web studies (Dennard et al. 2009;Ikeda et al. 2010;McClellan et al. 2010;Ward et al. 2010;Bond and Diamond 2011;Nosrati et al. 2014Nosrati et al. , 2018. Those models, for example, include SIAR using the Markov chain Monte Carlo method with an overall residual error term (Parnell et al. 2010;Parnell and Jackson 2013) and MixSIR (Moore and Semmens 2008), applying sample importance resampling. Recently, these models were applied to determine the quantitative contribution of different nitrate sources to nitrate contamination of groundwater and surface water (El Gaouzi et al. 2013;Korth et al. 2014;Zhang et al. 2015;Xu et al. 2016). In a hydrological context, similar models using the generalized likelihood uncertainty estimation (GLUE) methodology that also includes Markov Chain Monte Carlo methods, have been developed to better understand complex environmental systems (Beven and Freer 2001). Nevertheless, Bayesian models, such as SIAR, were rarely applied for describing microbial denitrification processes in published case studies to date (Yue et al. 2015;Xia et al. 2017;Li et al. 2019). Next to Bayesian models, Monte Carlos simulations are simpler and can also be used to model such processes coupled with statistical tools. Similar to Bayesian models, an advantage of Monte Carlo simulations is their inherent ability to characterize uncertainties and to provide probabilistic risk estimates of certain scenarios (Sadegh and Vrugt 2014). However, a best fit is highly dependent of the given data series and may also implicate uncertainties. To further advance isotope interpretation methods, Monte Carlo simulations can play an important role, especially for data from study sites with a complex hydrogeology and an input of different nitrate sources. Probability density functions can be assigned to each parameter reflecting uncertainty, and parameter sensitivity can be evaluated. Results can be evaluated in terms of probabilities, rather than deterministic values.
In a previous study (Wild et al. 2018), we analyzed stable isotopes of nitrate (yielding δ 15 N and δ 18 O) and concluded that nitrate contamination of groundwater near Hohenthann, Germany, may has originated from different nitrate sources comprising mineralization of soil organic nitrogen, as well as nitrification of ammonia and urea included in fertilizers. It was also stated by the authors that denitrification may have occurred in two of the wells screened in the perched, but not significantly in the main aquifer of the studied groundwater system.
To further assess and understand the processes influencing δ 15 N distributions in a groundwater system with a complex hydrogeological structure (main and perched aquifer structure), we ran inverse Monte Carlo simulations for the well-documented case mentioned above. In order to simulate δ 15 N-value distributions in groundwater arising from specific potential contributions, we proceeded stepwise by including (1) the land use and related input (agricultural versus non-agricultural land use and mixing of the nitrate sources manure, mineral fertilizer and precipitation), (2) hydrodynamic processes (advection and dispersion) in groundwater, and finally (3) possible microbial denitrification. Such contributions may explain δ 15 N-values observed in groundwater in more detail and may support the interpretation of isotope data, which have been analyzed in different simulation scenarios.

Methods
The study area is located near Hohenthann in Lower Bavaria, Germany, and it is affected by high nitrate concentrations in groundwater, reaching up to 120 mg/L (Wild et al. 2018). As displayed in Figure 1a, approximately 65% to 80% of the study area is agriculturally used with intensive hog farming and the cultivation of crops, predominantly maize. The subsurface is characterized by Quaternary and Tertiary clastic sediments dominated by sandy gravels, where a main aquifer is formed at around 45 m depth and deeper (Figure 1b which have been developed locally above clay layers. Hydrogeological characteristics and mean transit times (MTT) of groundwater, as identified in our previous study (Wild et al. 2018), suggest the definition of two hydrogeological units: the main aquifer (MA) and the perched aquifer (PA), where the latter combines all shallow perched aquifers investigated within the study area. Groundwater from the PA was sampled from shallow groundwater wells, in monthly sampling campaigns between December 2015 and March 2017. The MA was sampled only once (summer 2016) via deep groundwater wells (45 to 150 m bgl).
For both aquifer units, the PA and the MA, δ 15 N values of dissolved nitrate in groundwater were analyzed during the former survey (Wild et al. 2018  Note: Data for δ 15 N in manure, mineral fertilizer (MF) and precipitation (P) are taken from Kendall and McDonnell (1998), other parameters are defined for the aquifer system of the study area. PA, perched aquifer; MA, main aquifer; x, flow length; MTT, mean transit time of groundwater; ε, isotope enrichment factor; distr., distribution; min., minimum; max., maximum.
the observation period (December 2015 to March 2017) revealed to be low. Thus, considered δ 15 N values for both aquifers refer to the same time frame. Probability density functions were fitted to these measured δ 15 N values for the MA and PA, respectively, and δ 15 N frequency distributions determined, which were then compared to simulated isotopic distributions. The latter were generated by Monte Carlo simulations considering three scenarios that can be assumed for groundwater systems: (1) only mixing of different nitrate sources with characteristic δ 15 N signatures, (2) mixing combined with hydrodynamic processes (nitrate transport in groundwater without biotic or abiotic reactions), and (3) mixing and nitrate transport in groundwater affected by microbial denitrification. The three scenarios are described in the following, and Table 1 summarizes the considered parameter ranges and probability density functions. Scenario 1 simulates δ 15 N distributions for the mixing of different nitrate sources most relevant for the field site, including manure, mineral fertilizer and nitrate derived from precipitation. Typical δ 15 N distributions observed for these nitrate sources are reported by Kendall and McDonnell (1998), as shown in Figure S1 in the Supporting Information. These observed distributions were evaluated by the Anderson-Darling test (kurtosis sensitive), Kolmogorov-Smirnov test (sensitive to the mean of the distribution) and chi-squared (null hypothesis) test, and resulting p-values were compared in order to find best-fit distributions (obtained distributions are presented below) (Pettitt and Stephens 1977;Huber-Carol et al. 2008).
Agricultural (portion p 1 ) or non-agricultural land use (1-p 1 ) was considered, where the first can imply either the use of manure (p 2 ) or mineral fertilizer (1-p 2 ), and the latter is associated to nitrate derived from precipitation as the only nitrate source. δ 15 N values in groundwater (GW) are thus obtained as: Monte Carlo simulations applied random sampling of the fitted δ 15 N distributions for manure, mineral fertilizer and precipitation (δ 15 N M , δ 15 N MF and δ 15 N P ), where portions p 1 and p 2 were varied in steps of 0.05 (5%) for Scenario 1 and steps of 0.1 (10%) for further use in Scenario 2 and 3.
In Scenario 2 , it was assumed that nitrate released from the sources to groundwater is subject to hydrodynamic processes, while being transported along certain distances to the observation wells. As outlined in the introduction, a first modeling approach was done using analytical solutions that consider homogeneous conditions for the perched and the main aquifer. Otherwise multidimensional numerical advection-dispersion models could address a complex geological structure and hydrogeology, but such models need a detailed data set in high spatial resolution, which is not available for this site. Instead, in our study, aquifer heterogeneities were considered by random sampling from the probability density functions, which we have defined for the transport parameters of the analytical model (as described below). Our findings from modeling where then carefully compared to study site observations and literature findings. Accordingly an analytical solution for a 1D transport has been implemented based upon van Genuchten and Alves (1982), considering constant input to groundwater (at x = 0). δ 15 N values in groundwater as a function of time t, at location x downstream of the source, was modeled, accordingly, as: and u = v · (1 + 4kD/v 2 ) 0.5 , where δ 15 N GW is the initial δ 15 N in groundwater (input from sources at x = 0, Equation 1), v represents the groundwater flow velocity, x the flow length, t time and D the dispersion coefficient. The latter can be defined as D = α L · v x = P D · v · x , with longitudinal dispersivity α L and dispersion parameter P D . R is the retardation factor, and was set to 1 (no retardation assumed), and k is a first-order rate constant for degradation, set to zero in Scenario 2 (no degradation assumed). Probability density functions were defined for the Monte Carlos simulations, as described in the following. From 3 H/ 3 He measurements and modeling results, Wild et al. (2018) found a range of plausible P D values for the aquifer. To these values a normal distribution could be fitted with mean μ = 0.15 and standard deviation σ = 0.1, truncated by 0.01 and 0.3 (corresponding to minimum and maximum P D identified). The groundwater flow velocity v was calculated by dividing flow length x by the mean transit time (MTT). Corresponding to assumed ranges for the field site (average distance between nitrate sources and downstream groundwater wells), uniform distributions with x = 50 to 100 m and x = 500 to 1000 m were considered for the perched aquifer and the main aquifer, respectively. Based upon MTT determined from 3 H/ 3 He dating (Wild et al. 2018), an uniform distribution with 1 to 10 years was considered for the perched aquifer, and a log logistic distribution with location parameter (shift) γ = 4.18, scale parameter β = 18.01, shape parameter α = 2.60 was fitted for the main aquifer.
In Scenario 3 , hydrodynamic processes including microbial denitrification with isotopic enrichment of 15 N in the remaining nitrate were considered. For that, Equation 3 was used, where k (as part of coefficient u) was defined as k = μ · (α-1). There, μ [a −1 ] is a firstorder rate constant for microbial degradation of nitrate and α [−] is the isotope fractionation factor (derivation see Section S1 in the Supporting Information). The isotope enrichment factor ε [‰] is defined as ε = (α-1) · 1000, and a range of ε from −25 to −10‰ was considered, which has been observed for porous aquifers (Mariotti et al. 1981(Mariotti et al. , 1982Boettcher et al. 1990). This range was defined as a uniform distribution for the Monte Carlo simulations. For the first-order degradation rate constant μ, generic values of 0.1 and 1 a −1 were presumed in order to consider moderate and high microbial degradation in groundwater, exemplarily (based on typical ranges, Tesoriero and Puckett 2011).
Convergence was analyzed for all Monte Carlo simulations, where a relative stability of the calculated moments (average and variance) was reached after 5000 to 6000 trials, depending on the scenario and realization. This is shown qualitatively in Figure S32 in the Supporting Information for selected cases. We therefore decided to apply a slightly higher number of 10,000 trials for the Monte Carlo simulations, for which we applied the Microsoft Excel-Add-In @Risk (Palisade Decisiontools) as well as R version 3.5.1 (R Core Team 2018) implemented within RStudio 1.0.143 (RStudio, Inc.). Each realization of a scenario yielded random samples for the observation (using the fitted PDFs as described above) and the simulation.
Data pairs of 10,000 random samples, each, represented "the observation" and "the simulation" for every realization. We aimed at evaluating probabilities of δ 15 N values, that is, how observations could be explained by the processes considered in Scenarios 1 to 3. Thus, we set up histograms, in order to determine the frequency distribution of observed and simulated δ 15 N-values within certain ranges (bins). These bins were limited to −10 to +20‰ with an interval of 0.1‰. On the basis of these frequency distributions, the goodness of the model fit was evaluated by using the mean absolute error (MAE) and the coefficient of determination (R 2 ). The MAE indicates the mean absolute deviation between observation and modeling: where N is the number of bins (301 bins ranging from −10 to +20‰ with a constant bin width of 0.1‰), ζ ' i and ζ i is the observed and modeled frequency, that is, the number of random samples for observed and modeled δ 15 N, respectively, in each bin i . The smaller the MAE the better is the model fit. In this study, the MAE was preferred over the root mean squared error (RMSE), which is widely used in modeling studies, but often inappropriate and misinterpreted as it should only be applied for Gaussian distributions (Willmott and Matsuura 2005;Chai and Draxler 2014). As the best fit cannot be reduced to only the lowest MAE, we defined a best-fit range from the lowest MAE to the maximum acceptable MAE for each scenario. Each MAE relates to a specific run and thus to a specific realization of a scenario. The maximum acceptable MAE is defined as the fifth percentile of MAE (cumulative distribution of MAE for all realizations of a scenario).

Considered Nitrate Sources and Observed δ 15 N Values of Nitrate in Groundwater
Probability density functions (PDFs) were fitted to characteristic δ 15 N values of nitrate sources reported by Kendall and McDonnell (1998) (Figure S1 in the Supporting Information). A beta distribution was found as a best fit for manure (minimum = 3.25‰, maximum = 24.60‰, shape parameters α = 1.96 and β = 2.24) and normal distributions for precipitation and mineral fertilizer (with mean value μ = 0.62‰, standard deviation σ = 3.47‰ and μ = 2.06‰, σ = 2.00‰, respectively). These PDFs were used as input for the Monte Carlo simulations in order to define the characteristics of different nitrate sources. PDFs were also fitted to δ 15 N values observed in groundwater of the main aquifer and the perched aquifer, respectively, where logistic distributions could describe observations best (with location α = 6.199, scale β = 1.952 for the main aquifer and α = 9.221, β = 1.781 for the perched aquifer). Subsequently, Monte Carlo simulations were run (10,000

Scenario Best Fit range MAE [−] R 2 [−] M A E[ −] R 2 [−]
( Note: Values refer to realizations within the best fit range for all scenarios. Cf. Figure S27 for more information on MAE. trials) applying these PDFs in order to generate "measured" δ 15 N frequency distributions that could be compared to modeled δ 15 N values (as data pairs). Figure S2 in the Supporting Information shows measured values versus fitted distributions.

Scenario 1: Mixing of Different Nitrate Sources
Simulation results considering the mixing of possible nitrate sources reveal that a range of different portions concerning land use and related nitrate input could explain observed δ 15 N values in groundwater reasonably well. Those include agricultural (portion p 1 ) or non-agricultural land use (1-p 1 ), with manure (p 2 ) and mineral fertilizer (1-p 2 ) as nitrate sources for agricultural and precipitation for non-agricultural land use. Results are mainly discussed by means of cumulative frequency distributions and tile maps, but we also added some histograms in the SI to illustrate the resulting MAE and R 2 . A selection of good simulation curve fits is shown in Figure 2 for simulated cumulative frequency distributions of δ 15 N . "Good fits" were associated to a low range of MAE calculated for simulated versus observed δ 15 N distributions. This corresponds to MAEs from 8.47 to 10.16 for the perched aquifer and 5.27 to 9.67 for the main aquifer (Table 2, lowest MAE to 5th percentile MAE of all realizations for each aquifer, cf. Section 2 and Figure S27). Figure 3 presents tilemaps of R 2 and MAE for the considered realizations of Scenario 1, where blue to green colors indicate good fit, and red indicates bad fit.
For the perched aquifer, good curve fits were found with 70% to 100% of agricultural land use and 60% to 100% of manure application (Figures 2a and 3  was found for 100% agricultural land use and 60% manure (Table 2 and Table S1). However, such a high portion of agricultural land use is unrealistic: Burger (1993) estimated that 80% of the larger catchment area is agriculturally used, and from recent satellite images we estimated agricultural used areas to cover about 65% to 80%. Therefore, we can assume that 80% would be the maximum realistic percentage of agricultural land use, where the simulated realization with 75% agricultural land use and 85% manure can be seen as the best estimate (MAE of 8.49, with R 2 = 0.895). For the main aquifer, best fits were found with lower portions of agricultural land use between 40% and 60% (best estimate 45%) and relatively high portions of manure between 70% and 100% (best estimate 100%; Figures 2b and 3, right-hand side; Table 2 and Table S1).
The perched aquifer, located at shallow depths above 45 m bgl (meter below ground level), is characterized by relatively young groundwater with apparent MTT ranging from <4 years to 20 years. The deeper main aquifer extends from 45 to 150 m bgl, and it contains older groundwater with apparent MTT between 14 and 122 years (Wild et al. 2018). Results from the simulations of Scenario 1, considering the impact of possible nitrate source mixing on δ 15 N distribution, indicate higher portions of agricultural land use (75-80%) for the perched aquifer, compared to the main aquifer (40-60%, best estimate [b.e.] 45%). Manure seems to have contributed with 60% to 100% (b.e.: 85%) for the perched aquifer, and with 70% to 100% (b.e.: 100%) for the main aquifer. These findings point toward a change of land use within the past decennia, characterized by an increase of the agriculturally used area within the catchment. Although relative contribution of manure (usage of manure versus mineral fertilizer) seems to be constant or slightly lower than mineral fertilizer, the total amount of released manure seems to have increased with time. Thus, these simulation results might indicate a source-driven isotopic shift to heavier δ 15 N values of nitrate for the catchment area, away from less intensive farming with little livestock and mainly manure application (low use of mineral fertilizers) toward an increasingly intensive agricultural practice. This can be seen in Figure 2  be estimated at which proportion a specific source might have contributed to observed nitrate contamination in groundwater, as similarly done for other sites by applying Bayesian framework studies (El Gaouzi et al. 2013;Korth et al. 2014;Zhang et al. 2015;Xu et al. 2016).

Scenario 2: Hydrodynamic Processes
If hydrodynamic processes in groundwater (advection and dispersion) are considered, simulation results depend on the travel time of nitrate. As soon as the breakthrough of the isotopic signal, released at the source, has occurred at the observation well, good simulation curve fits were obtained ( Figures S3-S6, Table S2). Simulated δ 15 N frequency distributions (Figure 4 and Figure S30 in the Supporting Information, blue curves) are then similar to those obtained from Scenario (Sc.) 1 (Figure 2, as well as green curves in Figure 4 and Figure S30). Indeed, good simulation curve fits were obtained for the same range of source composition. For the perched aquifer, depending on the percentage of agricultural use and manure, simulations were within the best fit range (as defined in Table 2) after transport durations of 11 to 50 years ( Figures S3 and  S4, Table S2). Again, lowest MAE (6.37) was obtained for 100% agriculture and 60% manure (after 47 years), however this was not assumed realistic since not the whole catchment area is agriculturally used (see above). A more representative realization, considered as best estimate, was obtained for 80% agriculture with 80% manure after 50 years (MAE of 6.69). This source composition is similar to that yielding the best estimate for Sc. 1 (75% agriculture with 85% manure) and the MAEs of Sc.2 converge with time to the MAEs of Sc.1 as shown in Figures S15 and S16 in the Supporting Information.
For the main aquifer, best fit ranges were obtained after longer transport durations from 28 to 100 years ( Figures S5 and S6, Table S2). Associated source composition showed a wider range than for Sc. 1, with potentially 40% to 100% agriculture and 50% to 100% manure. The lowest MAE of 5.97 (R 2 = 0.977) was found for 50% agricultural land use with 100% manure after 65 years, so that, like for the perched aquifer, the source composition coincides relatively well with that of Sc. 1 (with 45% agriculture and 100% manure as best estimate). Thus, the time of breakthrough (transport duration) is an important unknown for Sc. 2, which needs careful consideration in order to derive realistic assumptions. It is mainly determined by advection (groundwater flow velocity v, in our case defined by the observed ranges and fitted PDFs for MTT and flow distance x), and it is also influenced by dispersion (PDF fitted for observed P D ). As soon as breakthrough has taken place, transport processes revealed only a low influence on the frequency distribution of δ 15 N. While comparing Scenario 1 and 2, the lowest MAEs are found for similar or even the same mixing portions in the perched and main aquifer ( Figures  S21 and S31 in the Supporting Information).

Scenario 3: Hydrodynamic Processes and Microbial Denitrification
In this scenario, mixing and transport (hydrodynamic processes) along with microbial denitrification were considered. First, we simulated a hypothetically high microbial activity, using a generic degradation rate constant μ of 1 a −1 . In comparison to transport without microbial denitrification, we can see a shift toward lower portions of agriculture and/or manure, which would allow similar simulated δ 15 N frequency distributions. This is illustrated in Figure 5, after a transport duration of 30 years for the perched aquifer and 60 years for the main aquifer. Best simulation results (low MAE) are shifted to   the left (less manure) and downwards (less agriculture) for Sc. 3, when comparing Figure 5a with 5c and 5b with 5d. The best fit range was found for 60% to 100% agriculture with 50% to 100% manure (after 10 years or longer) for the perched aquifer, and 30% to 100% agriculture with 20% to 100% manure (after 28 years or longer) for the main aquifer (Figures S7-S10, Table S3). This is a wider range for possible source compositions as compared to the previous scenarios. For the perched aquifer, similar combinations concerning the sources (portion of agricultural land associated with portion of manure), compared to Sc. 2, yielded good curve fits (Figure 4a, red curves). Lowest MAE (6.31) was again found for 100% agriculture and 60% manure after 16 years and thus earlier than for Sc. 2 (with 47 years) as displayed in Figures S17 and S18 in the Supporting Information. Again, since 100% agricultural land use is not representative for the study area, a portion of 80% agriculture with 80% manure and a transport duration of 46 years resulted in the most realistic estimate (MAE of 6.49, R 2 of 0.939). For Sc. 2 the best fit was obtained for 80% agriculture with 80% manure after 50 years, thus being very close.
The source composition was different for the main aquifer, where high portions of agriculture are associated with lower portions of manure, for obtaining similar δ 15 N frequency distributions (Figure 4b). The lowest MAE with 3.57 (R 2 = 0.984) was found for 60% agricultural land use with 40% manure after 81 years. Here we can see a clear difference compared to Sc. 2, where the best estimate was for a slightly lower portion of agricultural land use (50%), but much higher manure (100%), and after a shorter transport duration (65 years). The differences of the calculated MAE between Sc.2 and 3 are quite evident in Figure S23 in the Supporting Information.
The second assumption for Sc. 3, using a lower generic degradation rate constant μ of 0.1 a −1 , gave similar results compared to Sc. 2 (Figures S19, S20, S25, and S26 in the Supporting Information). In this case, microbial denitrification took only low influence on simulated δ 15 N frequency distributions (Figures S11-S14). For oxic groundwater, Tesoriero and Puckett (2011)  only be detectable for rates larger than 0.36 a −1 . This could explain the similarity to Sc. 2, which neglects microbial denitrification. For the perched aquifer, good estimates (as defined in Table 2) were obtained for 60% to 100% agriculture (slightly wider range as for Scenario 2) with 60% to 100% manure, after 11 to 50 years (as for Sc. 2). Lowest MAE was found for 100% agriculture and 60% manure (as for Sc. 2) after 32 years (for Sc. 2 it was 47 years). For the main aquifer, good estimates resulted from 40% to 100% agriculture (as for Sc. 2), however the manure best fit range started at a slightly lower percentage (40-100%). The best fit was found for 50% agriculture with 90% to 100% manure (similarly low MAE around 5.2) and thus very similar to Sc. 2. Simulation results show that microbial denitrification might have taken place in the perched aquifer, but it is rather unlikely for the main aquifer. For the perched aquifer, good simulation curve fits were obtained for Sc. 3 when considering transport and denitrification in groundwater combined with a nitrate input, which relates to expected portions of nitrate sources. This is also consistent with observations of Wild et al. (2018) that revealed significant denitrification in two shallow wells in the perched aquifer, but not in the residual shallow wells and springs. However, if assuming microbial denitrification for the main aquifer, good curve fits were predominantly found for unrealistic (unexpected) percentages of either agricultural use or manure. The main aquifer contains older groundwater (MTT of 14-122 years) compared to the perched aquifer (MTT <4-20 years, see above). Since the use of mineral fertilizers was lower in the past (Wild et al. 2018), we would expect a rather high percentage of manure. We can also assume a high percentage of agricultural use for the catchment area (80% was reported by Burger 1993). However, better curve fits prevailed for other source compositions (best estimate for 50% manure and 50% agriculture, where at least the latter is lower than expected). Therefore, the presence of microbial denitrification is less likely in the main aquifer as reasonable curve fits are not within a realistic range for the source composition (Table S3). These findings also agree with the calculated O 2 reduction rates and denitrification lag times of the investigated main aquifer, which suggest that it will take many decades to significantly reduce nitrate concentrations in the main aquifer via denitrification (Wild et al. 2018)

Conclusion
We investigated MC simulations as a decision support to interpret δ 15 N values of nitrate in groundwater. Different scenarios, such as mixing, combined with transport and microbial denitrification, were applied to study the influence of selected parameters on the evaluation of δ 15 N values in groundwater. Results show that the portion (percentage) of nitrate-releasing land use and specific nitrate sources in a catchment area along with the MTT of nitrate dissolved in groundwater are crucial factors when evaluating influences related to the mixing of different nitrate sources linked with transport and denitrification processes. However, if the sensitive parameters are well documented for a catchment area, MC simulations have the potential to support decision makers in the assessment of nitrate isotope data.