Multi-Tier Validation of a Macroscale Nitrogen Model for Groundwater Management in Watersheds Using Data from Different Monitoring Networks

Abstract

For the Hessian river basins, an area-differentiated modeling of the nitrogen input to the groundwater and surface waters was carried out for six diffuse input pathways and six point source input pathways on the basis of the geodata available at the state level. In this context, extensive plausibility checks of the model results were carried out using the data from several official monitoring networks at the state level. These include the comparison of modeled runoff components and input pathways for nitrogen using the data from the network of discharge monitoring stations. For the validation of the modeled nitrate concentrations in the leachate, the data from groundwater monitoring wells for controlling the chemical status of groundwater were used. The validation of the modeled nitrate inputs to the groundwater and denitrification in the groundwater was carried out using the data from a special monitoring network of groundwater monitoring wells that include N₂/Ar measurements. The data from the Surface Water Quality Monitoring Network were used to verify the plausibility of the modeled total N inputs to the surface waters from diffuse sources and from point sources. All of the model results evaluated by the plausibility checks prove that the nitrate pollution situation in Hesse is adequately represented by the model. This is a prerequisite for accepting the model results at the state level as a basis for developing and implementing regionally appropriate mitigation measures. The Hessian State Agency for Nature Conservation, Environment and Geology uses the model results in the broader context of the work on implementing the EU Water Framework Directive and the EU Nitrate Directive.

1. Introduction and Objective

Nitrogen models are used worldwide to simulate the nitrogen pollution of groundwater and surface waters at the regional scales [1,2,3,4]. Physically based reactive nitrate transport models, such as HYDRUS-1D [5] and the Daisy model [6], may be suitable for quantifying the nitrate loads and concentrations at the field scale [7], but their applicability is often limited by the regionally available database [8,9,10]. Risk assessment tools [11,12,13,14] require less-complex databases but are significantly limited in their conclusiveness, because nitrate pollution is often represented only in terms of relative degrees of contamination.

In order to be applicable for nitrogen management issues on the state scale, e.g., for assessing the dimension of the N reduction needed to achieve water quality targets, both the sources of nitrate pollution and the corresponding nitrate input pathways should be represented quantitatively in high spatial resolution [15]. Due to data limitations, risk assessment tools are applied in many parts of the world for this purpose [16,17,18,19]. In North America, coupled numerical modeling systems are applied to assess the compliance with quality standards and the impacts of agricultural activities on groundwater on the regional and national levels [17,20,21,22,23].

In Europe, the application of nitrogen flux models on the country level is often related to the requirements arising from the implementation of the EU Nitrate Directive [24], the EC Water Framework Directive [25] and the EU Marine Strategy Framework Directive [26], such as in Denmark [27], the Netherlands [28], France [29] and England [30].

In Germany, the development of nitrogen models at the federal state level started more than 20 years ago. Apart from the models Stoffbilanz [31] and MONERIS [32], the model system mGROWA-DENUZ-WEKU, which is used in 10 federal German states [33,34,35,36,37,38,39,40,41,42,43], has proven to be particularly successful in practice in Germany.

To ensure that the model results are included in political decision-making processes at the national level, it is essential that the model results are perceived by decision makers as reliable. An essential prerequisite for acceptance is the proof that the model results are in the same range of observed values at official monitoring networks at the state level. In the past two decades, there was an increase in the number of monitoring stations included in official monitoring networks used for measuring nitrate concentrations in Germany for their annual status reports to the European Environment Agency (EEA), for reporting under the EU Nitrates Directive (91/676/EEC) and for reviewing the chemical status of groundwater according to the EU Water Framework Directive (WFD 2000/60/EC). Accordingly, in recent years, the possibilities to check the plausibility of the model results determined with the model system mGROWA-DENUZ-WEKU using the time series data from such monitoring networks at the state level and, if necessary, to recalibrate the models, have increased [15,44,45,46].

The multi-tier validation presented here provides both an explicit evaluation of the performance of the individual modules of the mGROWA, DENUZ and WEKU model systems, and conclusions about the performance of the model systems with respect to the representation of the entire N flux, as follows:

• Do the runoff components simulated with the mGROWA model correspond to the discharge values from the Hessian network of gauging stations, so that it can be assumed that the modeled runoff adequately represents the regional discharge, and thus, the transport pathways for nitrogen?
• Is the nitrate concentration in the leachate modeled with the DENUZ model confirmed by the measured nitrate concentrations from the Hessian groundwater monitoring wells used for controlling the chemical status of groundwater, so that the related model result can be regarded as a reliable indicator for the regional nitrate pollution potential of groundwater?
• Is the nitrate degradation in the groundwater simulated in the WEKU model confirmed by the N₂/Ar measurements in the groundwater from a special monitoring network of groundwater monitoring wells that included the N₂/Ar measurements, confirming both the hydrochemical milieu characterization for the regional designation of denitrification conditions in the groundwater and the assumed denitrification kinetics in the groundwater?
• Do the modeled total N inputs to the surface waters correspond to the measured N loads at the Hessian Surface Water Quality Monitoring Network stations, confirming the overall performance of the model system regarding the modeled total N inputs to the surface waters from diffuse sources and point sources?

2. Methodology

The mGROWA-DENUZ-WEKU model system calculates the diffuse N inputs via the input pathways of erosion, wash-off, groundwater, natural interflow, drainage and atmospheric deposition to water bodies (Figure 1). Additionally, point source N inputs via municipal wastewater treatment plants, industrial/commercial direct discharges, from storm sewers in separate systems, through combined sewer overflows and from small wastewater treatment plants are assessed individually for each site and facility.

The modeling of diffuse N input was area-differentiated and carried out with high spatial resolution (100 m grid) for the entire Hesse state, so that every model calculation was performed for approx. 2.1 million grid cells. Dataof the mGROWA-DENUZ-WEKU model system are official state-wide and digitally available andof high quality in terms of actuality and consistency (Table 1).

The individual models of the mGROWA-DENUZ-WEKU model system, whose results are used in the multi-tier validation, were already described in detail in other publications, e.g., mGROWA [47,48], DENUZ [15,49] and WEKU [46,50], so only a brief summary of the models is given in this paper.

Table 1. Data type and data sources for the models mGROWA, DENUZ and WEKU.

Data Type	Data Source
Land use types	HLNUG 1: - ATKIS Basis-DLM25, - InVeKoS 2011–2016
Agricultural data	Thünen-Institute: - N balance surpluses from RAUMIS model [51]
Atmospheric N deposition	German Environment Agency (UBA): - N deposition from PINETI-3 model
Imperviousness	Copernicus Land Monitoring Service: - High Resolution Layer Imperviousness (2015) 20 m, https://land.copernicus.eu/pan-european/high-resolution-layers/imperviousness (accessed on 6 June 2023)
Digital elevation model	Federal Agency for Cartography and Geodesy (BKG): - Digital elevation model, resolution 25 m (DGM25)
River system	Federal Agency for Cartography and Geodesy (BKG): - Digital Landscape Model 1:250.000 (DLM250)
Soil	HLNUG 1: - Soil map BFD 1:50.000 - Soil horizon parameters: thickness, field capacity, effective field capacity, bulk density, fine soil type, fine soil group, substance volume - Soil profile parameters: water logging tendency, pedological depth to water table
Drainage areas	Newly derived - Based on methodology in [52]
Erosion data	HLNUG 1: - ‘ABAG’—Factors from Erosion Atlas 2018
Climate data	Climate Data Center (CDC) of the German Weather Service (DWD): - Precipitation and potential evapotranspiration for hydrological years 1991–2020
Hydrogeology	From various studies: - Hydrogeological areas [53] - Base flow indices [54] - Hydraulic conductivity [53] - Groundwater surface [53]
Contents of total phosphorous in topsoil	From study [55]
Groundwater quality data	HLNUG 1: - Fe, Mn, NO3, O2, DOC, (2010–2020)
Runoff and river water quality data	HLNUG 1: - Catchment areas of gauges and water quality measuring sites - Runoff time series (1991–2020) - N time series

Note(s): ¹ Hessian Agency for Nature Conservation, Environment and Geology.

Table 1. Data type and data sources for the models mGROWA, DENUZ and WEKU.

Data Type	Data Source
Land use types	HLNUG 1: - ATKIS Basis-DLM25, - InVeKoS 2011–2016
Agricultural data	Thünen-Institute: - N balance surpluses from RAUMIS model [51]
Atmospheric N deposition	German Environment Agency (UBA): - N deposition from PINETI-3 model
Imperviousness	Copernicus Land Monitoring Service: - High Resolution Layer Imperviousness (2015) 20 m, https://land.copernicus.eu/pan-european/high-resolution-layers/imperviousness (accessed on 6 June 2023)
Digital elevation model	Federal Agency for Cartography and Geodesy (BKG): - Digital elevation model, resolution 25 m (DGM25)
River system	Federal Agency for Cartography and Geodesy (BKG): - Digital Landscape Model 1:250.000 (DLM250)
Soil	HLNUG 1: - Soil map BFD 1:50.000 - Soil horizon parameters: thickness, field capacity, effective field capacity, bulk density, fine soil type, fine soil group, substance volume - Soil profile parameters: water logging tendency, pedological depth to water table
Drainage areas	Newly derived - Based on methodology in [52]
Erosion data	HLNUG 1: - ‘ABAG’—Factors from Erosion Atlas 2018
Climate data	Climate Data Center (CDC) of the German Weather Service (DWD): - Precipitation and potential evapotranspiration for hydrological years 1991–2020
Hydrogeology	From various studies: - Hydrogeological areas [53] - Base flow indices [54] - Hydraulic conductivity [53] - Groundwater surface [53]
Contents of total phosphorous in topsoil	From study [55]
Groundwater quality data	HLNUG 1: - Fe, Mn, NO3, O2, DOC, (2010–2020)
Runoff and river water quality data	HLNUG 1: - Catchment areas of gauges and water quality measuring sites - Runoff time series (1991–2020) - N time series

Note(s): ¹ Hessian Agency for Nature Conservation, Environment and Geology.

2.1. The mGROWA Model

To determine the amount of nitrate entering surface waters via runoff and erosion processes, the expected long-term mean surface runoff is first calculated as a fraction of the total runoff [56]. Based on the mGROWA model [47,48], which is a deterministic, conceptual, grid-based, area-differentiated water balance model, first the actual evapotranspiration and subsequently the direct runoff components, interflow and the drainage runoff, and runoff from urban areas as well as the groundwater recharge, are simulated. In this way, the spatial variability of hydrological conditions in Hesse and the associated regional differences in the relevance of the individual runoff components for the N input are represented.

The basis of the mGROWA model is the water balance equation with its climate, runoff and storage variables. To represent a long-term mean—and thus, a regionally typical hydrological situation—the water balance modeling was performed for a hydrological period of 30 years, more precisely, the period of 1991–2020. The modeling was performed in two steps (Figure 2).

In the first step, the actual evapotranspiration and runoff generation were calculated in daily time steps. Actual evapotranspiration was calculated based on potential evaporation above pasture [57,58], land-use-specific coefficients, a topographic correction function [59] and the Disse equation to account for the dependence of evapotranspiration on soil moisture [60]. To calculate soil moisture and leachate dynamics for sites covered with vegetation, the multi-layer soil water balance model BOWAB [61] was integrated into mGROWA. Depending on the water tension in the root zone, a capillary rise was calculated for groundwater-affected soils, i.e., at such sites, evapotranspiration from groundwater was determined. Additionally, mGROWA calculated the water balance for urban areas [62].

In the second step, the separation of the leachate into groundwater recharge, and the direct runoff components urban direct runoff, interflow and drainage runoff were carried out in monthly time steps. For this purpose, the site characteristics relevant for the separation of the leachate were identified and parameterized. In areas with unconsolidated rock, the installation depth of agricultural drainage systems [38,52] and the seasonally fluctuating groundwater tables control the separation of the leachate into drainage runoff and groundwater recharge. This means that drainage runoff is formed exclusively when the groundwater table has risen above the installation depth of the drainage systems.

For bedrock regions, the leachate rate was separated into interflow and groundwater recharge using calibrated BFI values [63] that reflect the ratio of groundwater recharge to total runoff, which is primarily differentiated by the hydraulic conductivities of the groundwater bearing rocks. Interflow is considered as the direct subsurface runoff component that does not contribute to groundwater recharge [64,65]. The BFI concept assumes that groundwater discharge and groundwater recharge of a basin are equal on a long-term average, provided that the total groundwater resources of the area remain unchanged. BFI values were already determined in Europe for a number of characteristic site properties [59,66,67,68,69] and are also available for Hesse on a grid basis [54].

2.2. The DENUZ Model

Determination of diffuse N outputs from soils and the nitrate concentration in the leachate is simulated with the reactive N transport model DENUZ [49]. DENUZ considers agricultural N balance surpluses and the atmospheric N deposition as diffuse nitrogen sources (Figure 3).

One part of the diffuse N input is potentially immobilized in the soil or incorporated into the biomass. N immobilization is assumed to be temperature dependent, according to Nagel and Gadow [70], for soils of semi-natural areas, i.e., mainly forest areas and non-agricultural grassland and heathland. For agriculturally used grassland, the amount of N immobilized in the soil is assumed as a percentage of the total N input to the soil [71].

For the part of the nitrogen stored in the biomass of forests (“uptake into the stand of forest areas”), fixed values are used for deciduous, coniferous and mixed forests [72,73]. For arable soils, it is assumed that the soils are N saturated due to years of fertilization and tillage [74], so that N contents in the topsoil remain nearly constant. The N outputs from arable soils then correspond to the amounts of N not taken up by the crop (mainly N surpluses from the fertilization and atmospheric N deposition) minus denitrification losses in the root zone of the soil.

Denitrification in soil is a process that occurs in all soils; however, the extent and kinetics of denitrification in soil depend on many different influencing factors. According to reviews on the literature [75,76], denitrification losses in soil occur mainly in the root zone at low oxygen and large water contents, as well as large contents of organic matter. In contrast, low denitrification rates are expected in well-aerated soils with low residence time of leachate in soil and low water and organic matter contents. In the DENUZ model, denitrification in soil and the corresponding N output from soil after residence time (N(t)) are simulated based on Michaelis–Menten kinetics as follows:

$$\frac{dN\left(t\right)}{dt}+D_{max}·\frac{N\left(t\right)}{N\left(t\right)+k}=0$$

N(t)	N output from soil after residence time t	(kg N/(ha·a))
t	Residence time of leachate in soil	(a)
Dmax	Maximum denitrification rate	(kg N/(ha·a))
k	Michaelis constant	(kg N/(ha·a))

The maximum denitrification rate in the root zone of soils (D_max) depends on the influence of groundwater and waterlogging, the initial geological substrate and the soil type [75,77]. The residence time of leachate in the soil (t) has a significant influence on the denitrification rate [78,79,80]. To account for this in the DENUZ model, D_max-values are related to the residence time of leachate in the soil. This ensures that denitrification is not overestimated in sandy soils displaying short residence times and underestimated in loamy soils displaying long residence times. The residence time of the leachate in the soil is calculated using the integrated effective field capacity of the soil in the root zone and the leachate rate [81] as follows:

$$t=\frac{nFK·We}{q_{Sw}}$$

t	Residence time of leachate in soil	(a)
qsw	Leachate rate	(mm/a)
nFK	Effective field capacity	(mm/dm)
We	Effective rooting depth	(dm)

The N output from soil reaches surface waters coupled to the runoff components determined by the mGROWA model (see above). One part reaches the surface waters coupled to the direct runoff components (surface runoff, drainage runoff, natural interflow, urban direct runoff) without consideration of further denitrification processes, and another part reaches the aquifer via groundwater recharge.

2.3. The WEKU Model

Nitrogen output from soil that does not reach surface waters via direct runoff enters the aquifer via groundwater recharge. In the WEKU model [82,83], the N input into the groundwater is assumed to be a proportion of the N output from soil equal to the proportion of groundwater recharge in the total runoff. In addition to the soil-borne N inputs to groundwater, N emissions from leaking urban systems and small wastewater treatment plants discharging into groundwater are considered as additional N sources to groundwater.

In groundwater, nitrate may be completely [50] or partially [8] denitrified on its way to a receiving surface water or to a groundwater monitoring well. The prerequisite for denitrification in the aquifer is reducing milieu conditions. If these are present, denitrification in the aquifer is determined in the WEKU model using first-order reaction kinetics [84,85] depending on the nitrate input into the aquifer, travel times and the nitrate degradation conditions in groundwater [86,87]. This is in accordance with the reaction kinetics determined for sites in the Netherlands based on extensive field investigations over several years [79,80,81,82,83,84,85,86,87,88,89].

$$\frac{dN\left(t\right)}{dt}+k_{n}N(t)=0$$

N(t)	Nitrate content in groundwater after travel time in aquifer	(kg N/(ha·a))
t	Travel time of groundwater	(a)
kn	Denitrification constant	(a−1)

The modeling of travel times in groundwater comprises three steps. In the first step, the hydraulic gradient and the receptor of groundwater discharge (surface waters, monitoring wells) are derived from digital elevation models of the groundwater surface. In the second step, groundwater velocity is calculated, which characterizes the movement of a water particle between two points in the flow direction of groundwater. The parameters needed to calculate groundwater velocity are hydraulic conductivity, the effective void ratio and the hydraulic gradient. While the values for hydraulic conductivity were taken from a hydrogeological map, the effective void ratios were estimated. For unconsolidated rocks, a correlation of hydraulic conductivity values and effective void ratios was used [90], while for bedrock regions, Hessian-specific values were applied [53]. In the third step, the travel time for the individual grids is calculated from the groundwater velocities determined for the entire flow path into a receiving water or groundwater monitoring well.

Denitrification conditions in the aquifer are determined by a hydrogeochemical milieu characterization based on groundwater quality data [46]. The concentrations of the redox-sensitive parameters (iron, manganese, nitrate, oxygen, DOC) were first interpolated using the IDW method. These interpolated concentrations were then assigned point values, with the highest class of the respective contents lying in a range that is typical for reduced aquifers, and thus, nitrate-degrading aquifers. In the lowest class, on the other hand, the respective contents are in a range that is typical for oxidized aquifers, and thus, non-denitrifying aquifers.

To provide an indication of the nitrate degradation capacity in groundwater from the regional interaction of the five individual redox-sensitive parameters, the point values of the regionalized and classified individual parameters were summed for each grid cell (100 × 100 m). These point values were then assigned a denitrification capacity, i.e., half-lives or reaction constants of denitrification, so that they could be integrated into the WEKU model for quantifying denitrification in the aquifer and groundwater-borne nitrate inputs to surface waters and monitoring wells, respectively. For this purpose, reaction constants/half-lives from the literature were used [33,40,50,85,89,91,92,93,94,95,96,97,98,99].

The models described in the Section 2.1, Section 2.2 and Section 2.3 are each implemented in standalone software written in Java. Graphical representations in this study were generated using RStudio 2023.03.1 (Figures 4, 5 and 8), ArcMap 10.8.2 (Figure 6), Java (Figure 7) and MS Excel version 2304 (Figures 9 and 10).

3. Validation of Modeled Runoff Components

The reliability of the modeled runoff determined by the mGROWA model is checked against the observed discharge values (MQ) from the gauging stations. For this purpose, the area-differentiated runoff rates modeled in mGROWA in 100 m × 100 m grids are integrated over the hydrological catchments of the gauging stations.

After excluding from the Hessian network of discharge gauging stations the gauging stations for which observation data was available for less than 20 years, and whose catchment areas were located to less than 85% in the Hessian state territory, in the end, the discharge data of 78 gauging stations were suitable for the plausibility check.

The model performance was evaluated using the quantitative statistical parameters NSE (Nash–Sutcliffe efficiency) and PBIAS (percent bias). The NSE [100] indicates the approximation of the simulated and modeled discharges to the 1:1 line. The NSE values range from zero to one. The closer the NSE value is to one, the better the model performance. The PBIAS value [101] indicates the tendency of the simulated runoff to be underestimated or overestimated. The closer the value is to zero, the less the model tends to overestimate or underestimate the measured total runoff. The calculation of the NSE and PBIAS was made on a weighted basis, whereas the larger catchment areas were weighted more heavily [102]. Figure 4 shows the modeled total runoff and the agreement at the gauging stations.

As indicated by an NSE of 0.88 and a PBIAS of −2.9%, the adjustment of the modeled total runoff rates to the observed MQ values at the gauging stations represents a very good agreement. In the context of the modeled nitrate inputs to the groundwater or surface waters, it is crucial whether the model adequately represents the input pathways for nutrients, i.e., the runoff components. The runoff components depicted by the mGROWA model include, on the one hand, the direct runoff components (direct runoff from agricultural drainage systems, interflow, direct runoff from settlement areas and surface runoff) and, on the other hand, baseflow, which largely corresponds to the groundwater recharge on a long-term average.

The possibility of validating the modeled direct runoff rates at the state level is limited by the lack of official monitoring networks for the long-term recording of direct runoff components, such as the drainage runoff. The only runoff component for which plausibility checks are usually performed is the groundwater recharge. For areas in Germany where unconsolidated rocks dominate, the mean groundwater runoff is often approximated by the minimum monthly low flow (MoMNQ) [103]. For Hesse, this applies to only one gauging station. All other Hessian catchments are dominated by consolidated rocks, for which the method proposed by Demuth in [104] is suitable for runoff separation. It uses a permanent curve to determine the mean baseflow, yielding either an S-shaped curve type with a linear range, or a parabolic curve type. Since the Demuth method did not provide reasonable results for 23 Hessian catchments with a parabolic curve type, these results were not used for the evaluation. Figure 5 shows the groundwater recharge rates modeled with the mGROWA model (left) and the agreement with the observed baseflow at the remaining 55 river gauges (right).

The plausibility check of the modeled groundwater recharge rates shows an NSE of 0.26 and a PBIAS of 9.6%. The latter indicates a slight underestimation of the modeled values, but overall, an acceptable model performance was achieved. When evaluating the results of the plausibility check of the groundwater recharge, it should be noted that the scattering effects of the three methods (total runoff calculation, groundwater recharge calculation, runoff separation from discharge) overlap. Therefore, the adjustment of the modeled groundwater recharge rates to the separated baseflow rates is always somewhat lower than the adjustment of the modeled total runoff rates to the measured MQ values.

The plausibility checks of the mGROWA model results in numerous Hessian catchments, which allows for several conclusions to be drawn about the subsequent modeling of the nitrate inputs. The only runoff component that represents an important input pathway for nitrate and that was explicitly validated in the plausibility checks is the groundwater recharge. The discharge data from the river gauges are of limited use to validate individual direct runoff components because all direct runoff components enter surface waters relatively soon after a precipitation event and are difficult to distinguish from each other. Still, the sum of the direct runoff components is equal to the difference between the total runoff and groundwater recharge, so in this way, the validity of the sum of all modeled direct runoff components was also confirmed.

Since the two main components of direct runoff, i.e., interflow and drainage runoff, are largely mutually exclusive in a region, this can be regarded as a validation of the modeled interflow and drainage rates in that region. Thus, the validation has shown that the major runoff components of the groundwater recharge, interflow and drainage runoff can be assumed to be regionally representative. Accordingly, it can be concluded that the main input pathways for nitrate were modeled in a regionally realistic and representative manner.

The same applies to the modeled leachate rates leaving the root zone of the soil vertically. The most reliable method for validating the modeled leachate rates is to compare them to the leachate rates from the lysimeter measurements [105]. However, this type of validation is not practical for a state-wide assessment because there are generally not enough lysimeter stations available to provide statistically significant results for the site combinations that occur in a state.

In the hydrological water balance, the leachate rate corresponds to the amount of water that is neither lost due to evapotranspiration nor discharged as surface runoff. Since the share of surface runoff in the total runoff in Germany is insignificant at the catchment level [106], the leachate rates are generally of the same order of magnitude as the total runoff rates. Accordingly, the plausibility check of the total runoff (Figure 4) indirectly validated the modeled leachate rate. This fact is significant because the leachate rate is an important parameter for calculating the nitrate concentrations in the leachate and assessing the nitrate pollution risk in a region.

4. Validation of Modeled Nitrate Concentration in the Leachate

Authors should discuss the results and how they can be interpreted from the perspective of previous studies and of the working hypotheses. The findings and their implications should be discussed in the broadest context possible. Future research directions may also be highlighted.

Since there is no obligation in Germany to document and record the use of fertilizers at the farm level [107], there are hardly any records of the locally applied amount of N fertilizer and the resulting N surpluses, which would facilitate the verification of the N balance surpluses modeled by the RAUMIS model [108,109]. Indirectly, however, such a plausibility check is possible, precisely by means of the modeled nitrate concentration in the leachate, determined on the one hand by the displaceable N load in the soil originating from agricultural N balance surpluses and atmospheric N deposition, and on the other hand, by the leachate rate diluting this N load. As denitrification processes in the unsaturated zone below the soil zone can be neglected, the nitrate concentration in the leachate can be regarded as the key indicator of the nitrate pollution of groundwater [15]. Additionally, it is an indication of the nitrate concentrations reaching surface waters via the direct runoff components, interflow and drainage runoff, respectively.

An extensive plausibility check of the magnitude and spatial distribution of the modeled nitrate concentration in the leachate is therefore indispensable. The measured values from the soil depth profiles, suction probes and lysimeters are very well suitable for checking the plausibility of the modeled nitrate concentrations in the leachate. However, the number of available monitoring stations is not high enough to draw conclusions about the validity of the model results at the state level [46].

For this reason, a preselection method is applied to identify suitable groundwater monitoring wells from the upper aquifer and springs whose observed nitrate concentrations may be compared to the modeled nitrate concentrations in the leachate [46]. The preselection criteria are the exclusive inclusion of monitoring wells filtered near the surface (up to 30 m below GWO), as well as the exclusive inclusion of groundwater monitoring wells and springs indicating an oxidative milieu. By applying these preselection criteria to the data from the network of the Hessian groundwater monitoring wells and springs for controlling the chemical status of groundwater, 2256 groundwater monitoring wells and springs were identified as suitable for the plausibility check. For these groundwater monitoring wells and springs, the mean values are derived from the individual nitrate measurements and compared to the modeled mean nitrate concentration in the leachate of the inflow area of the corresponding measuring point or spring.

Figure 6 shows an overview of the nitrate concentrations in the leachate modeled across Hesse in the 100 m × 100 m grids and the preselected monitoring wells and springs with nitrate concentrations in groundwater (dots). Since the same class widths and the same color gradation were chosen, the values are visually directly comparable. It becomes apparent that the spatial patterns are reflected very well, i.e., the regions with high nitrate concentrations in the groundwater and springs are just as well represented by the modeled nitrate concentrations in the leachate as the regions with low nitrate concentrations in the groundwater and springs.

To achieve a systematic, comprehensible assessment of model validity, the difference between the absolute values of the modeled nitrate concentrations in the leachate and the observed nitrate concentrations in the groundwater and springs was determined and subsequently evaluated statistically. Ideally, all the differences should be zero; in practice, however, the differences are distributed around a mean value. The closer the difference between the modeled and observed values is to “0”, the better the agreement in general, and the smaller the distribution is, the better the correlation of the measured and modeled values is in detail.

The frequency distribution over all 2256 monitoring stations in Figure 7a shows that half of the modeled nitrate concentrations in the leachate deviate less than 0.17 mg NO₃/L from the observed values, and 75% less than 7 mg NO₃/L. This can be considered as a good agreement, which is confirmed by the standard deviation of 21.5 mg NO₃/L. The scatter shows a relative normal distribution, indicating that the modeled nitrate concentrations in the leachate show neither a tendency to overestimate nor to underestimate.

When evaluated separately according to the three main types of land use, namely, arable land, grassland and forest, a more differentiated result emerges. Especially for the land use categories forest (Figure 7d) and pasture (Figure 7c), displaying 1196 and 705 monitoring stations, respectively, the difference between the simulated nitrate concentrations in the leachate and the observed nitrate concentrations was quite low, as indicated for the forest sites by a P50 of 0.48 mg/L, a P75 of 3.87 mg NO₃/L and a standard deviation of 15.65 mg NO₃/L, and for the grassland sites by a P50 of 0.94 mg NO₃/L, a P75 of 8.3 mg NO₃/L and a standard deviation of 16.23 mg NO₃/L.

As indicated by a P50 of 7.22 mg NO₃/L, a P75 of 21.78 mg NO₃/L and a standard deviation of 42.52 mg NO₃/L, the difference between the simulated nitrate concentrations in the leachate and the observed nitrate for the land use category arable land (Figure 7b) is significantly higher. This scattering of values is probably due to small-scale differences in the cultivated crops or in the farmer-specific level of fertilization, which is presently not considered in the calculation of agricultural N balance surpluses with the RAUMIS model on the community level. This is supported by the model calculations of the nitrate concentrations in the leachate for the regions in which the deviation could significantly be reduced in case the generalized N balance surpluses at the community level were replaced by local (plot-based) N balance surpluses [40]. Against this background, the high scatter of values is currently no argument for recalibrating the DENUZ model with the goal of bringing the modeled values closer to the measured values. Instead, it underpins the need to determine the agricultural N balance surpluses with a higher spatial resolution and accuracy.

The overall good correspondence of the spatial patterns of the modeled nitrate concentrations in the leachate with the observed nitrate concentrations in the groundwater and the low deviation of values for the forest and grassland sites has several implications. Firstly, it is proof that the parameterization of the DENUZ model for accounting the N fixation in the soil via immobilization and uptake into the stand as well as for accounting denitrification in the soil was performed correctly.

Secondly, this is proof that it is safe to use the modeled nitrate concentration in the leachate as a reference for regional nitrogen management issues. This includes the designation of hotspot areas of groundwater contamination and the assessment of the regional N reduction requirement to ensure a nitrate concentration in leachate of less than 50 mg NO₃/L [34]. Therefore, it is justified to use the identified N reduction requirement as a reference to simulate the effects of specific N reduction measures on groundwater nitrate contamination [110].

5. Validation of Modeled Denitrification Rates in Groundwater

The denitrification in the groundwater modeled with the WEKU model was compared with the denitrification quantified at the groundwater monitoring wells. The latter was based on the N₂/Ar measurements at the groundwater monitoring wells, from which the N₂ excess in the groundwater was determined [111,112], to assess from this the extent of nitrate degradation in the groundwater and the initial nitrate concentration of the groundwater recharge. The N₂ excess is calculated from the measured concentrations of molecular nitrogen (N₂) and argon (Ar) in the groundwater samples [113]. In contrast to molecular nitrogen, argon, as a noble gas, is not submitted to biochemical transformation processes in groundwater. Given a temperature of recharged groundwater of about 10 °C, which is typical for Germany [114], and an atmospheric composition like that of air in the unsaturated zone, the concentration of dissolved molecular nitrogen in the groundwater is 17.7 mg N₂/L and 0.67 mg Ar/L in the case of argon [115]. The denitrification in the aquifer increases the concentration of molecular nitrogen, shifting the measured N₂/Ar ratio to higher values. The excess N calculated in this way can then be attributed to the nitrate denitrified in groundwater.

For the state of Hesse, the N₂/Ar measurements were available for 278 monitoring sites from a special monitoring network of groundwater monitoring wells that included N₂/Ar measurements for the period of 2018–2022. After excluding the N₂/Ar monitoring wells that were filtered deeper than 30 m below the groundwater surface, and N₂/Ar monitoring wells for which no inflow areas could be identified, 151 N₂/Ar monitoring wells remained for the plausibility check. For each of these wells, the denitrification was determined from the N₂/Ar ratio and then compared to the mean denitrification in the groundwater calculated by the WEKU model for the inflow areas of the monitoring wells.

To compare the modeled denitrification with the denitrification calculated from the N₂/Ar measurements, four classes were formed (0–25%, 25–50%, 50–75%, 75–100% denitrification of initial nitrate in groundwater) (Figure 8, left). According to this classification, there is a good agreement for 54% of the sampling stations, i.e., simulated and measured nitrate degradation fall into the same class, and another 26% show a satisfactory agreement, i.e., a deviation by one class (Figure 8, top right). Most of the monitoring wells that show a good agreement with the modeled values are in the 75–100% and 0–25% denitrification classes. Thus, the regions where the denitrification in the groundwater is low appear to be represented just as well in the model as the regions where the denitrification is high. In contrast, the pairs of values with medium-to-high denitrification (25–75%) are in the minority. This is an indication that denitrification in groundwater is a process that either occurs entirely or not at all.

Figure 8 (bottom right) shows that the model results tend to slightly underestimate the measured denitrification in the groundwater. However, since there is no discernible spatial clustering of the N₂/Ar stations with good or poor agreement, the discrepancies likely reflect small-scale heterogeneities in the WEKU model input data. These include, for example, regional blurring with respect to the N balance surpluses calculated at the community level, and with respect to the groundwater monitoring well inflow areas derived from the digital groundwater surface available on the state level.

Based on this result, we conclude that both the reliability of the derived denitrification conditions and the associated reaction constants of denitrification in the groundwater are confirmed. With regard to the modeled nitrate input from the groundwater to the surface waters or monitoring wells, the results confirm the adequate simulation of this important input pathway.

6. Validation of Modeled Total N Inputs to Surface Waters

In the previous sections, the results of the plausibility checks for the individual modules of the mGROWA-DENUZ-WEKU model system were presented. This section deals with the result of the plausibility check of the modeled total N inputs to the surface waters from diffuse and point sources. Figure 9 provides an overview of the contribution of the individual diffuse and point N inputs to the total N input to Hessian surface waters in t/a, as well as the respective percentage of the individual N input pathways to the total N inputs to the receiving waters.

Figure 9 shows that the N input into the receiving waters mainly occurs via the diffuse input pathways “interflow” and “groundwater”. About 67% of the total nitrogen input (32,894 t N/a) enters the receiving waters via these two input pathways. All other diffuse N input pathways (drainage, erosion, runoff, deposition on water surfaces) are of minor importance for the N load of the Hessian receiving waters. The point source N inputs are dominated by the municipal wastewater treatment plants (approx. 7187 t N/a). With a ratio of point to diffuse sources of 30% to 70%, however, the N input from diffuse sources represents by far the major part.

Apart from the plausibility check of groundwater-borne N inputs based on the data from the N₂/Ar measurements, there is a lack of monitoring data in Hesse at the state level to validate the modeled N inputs for the individual input pathways. Accordingly, only the sum of the modeled N inputs to the surface waters could be verified using the measured N loads in surface waters. For this purpose, the long-term monitoring data (discharge and N concentration) of 98 receiving waters of the Hessian Water Quality Monitoring Network were available.

For this purpose, the N inputs to the surface waters from diffuse sources and from point sources were aggregated for the 98 catchments concerned. In addition, the N retention in the surface waters was considered based on the algorithms published by Behrendt and Opitz [116]. For the 98 receiving waters of the Hessian Water Quality Monitoring Network, the simulated N retention in the surface waters (between 16% and 60% of the total N input) was subtracted from the sum of the modeled N inputs. The remaining N quantity is assumed to be identical to the N load in the receiving waters and was accordingly compared with the observed N loads of the 98 receiving waters of the Hessian Water Quality Monitoring Network.

For these 98 gauge-related catchments, the plausibility check resulted in a coefficient of determination of R² = 0.99, as well as an NSE of 0.98 and a PBIAS of −4.0% (Figure 10). This represents a good agreement in the context of the available data and the size and heterogeneity of site conditions in the federal state of Hesse and can be regarded as a confirmation of the overall good performance of the mGROWA-DENUZ-WEKU model system for determining the N input to the surface waters from diffuse and point sources.

The analysis of N inputs by input pathways shows that the two most important input pathways in Hesse, interflow and groundwater, dominate in many rural catchments. In catchments with a higher population density, the N contribution from point sources increases. The overall good agreement at the monitoring sites in rural areas as well as in urban areas indicates that the interplay of the individual N input pathways was modeled realistically and representatively.

7. Discussion of Multi-Tier Validation

Checking the plausibility of the model results against the measured values is part of every modeling process. In the case of the N model system mGROWA-DENUZ-WEKU used here to support the N management at the state level in Hesse, the aim is to demonstrate that the current nitrate pollution of the groundwater and surface waters in the various Hessian regions is adequately represented. For this purpose, all results of the individual models, for which the measurement data with the corresponding measurement frequency were available state-wide, were subjected to a plausibility check.

The total runoff rates and groundwater recharge rates modeled with mGROWA were validated using the runoff data observed at the gauging stations. For this purpose, the values modeled on 100 m × 100 m grids for the period of 1991–2020 are summed for gauging the station-related catchments and compared with the observed runoff records of this period at the respective gauging stations. It should be noted, however, that a good agreement, indicated by an NSE or an R2 value close to one, does not prove that the values modeled for each 100 m × 100 m grid are correct.

Theoretically, it is possible that the over- and underestimations in the grid-wise modeled runoff components balance each other out and pretend to have a good agreement in the NSE and R2, respectively. If—as in the case of the modeling carried out here—the goal is to represent the runoff components as spatially high-resolution input pathways for nitrate, the corresponding plausibility check poses a fundamental problem, since there are no monitoring systems at the state level for recording and balancing individual runoff components at discrete locations.

Nevertheless, to assess the representativeness of the modeled runoff components/input pathways, it is important to include the runoff records of as many gauge-related catchments as possible, and to determine the PBIAS in the model validation. This is because the more catchments are included in the comparison and the more of these catchments have a PBIAS value close to zero, the more likely it can be assumed that the model results for the individual grids are not due to a compensating over- or underestimation of the modeled runoff values.

Against this background, the 79 gauge-related catchments that were available for the plausibility check of the total discharge levels, as well as the 55 gauge-related catchments that were available for the plausibility check of the groundwater recharge levels, are a solid basis for the assessment of the model performance. Since the runoff components determined by the mGROWA model show a good agreement at the gauging stations, it is safe to assume that the regional N emission pathways are sufficiently well represented in the model.

The modeled nitrate concentrations in the leachate were confirmed at 2256 groundwater monitoring wells from the upper aquifer and springs. However, it should be remembered that the modeled nitrate concentrations in the leachate and the measured nitrate concentrations in the groundwater have a different reference period. While the modeled values represent long-term averages, the measured values reflect the situation at a certain time (random sample). In addition, there are different spatial reference levels here, namely, the soil below the root zone for the modeled nitrate concentrations in the leachate, and the shallow groundwater for the nitrate concentrations in the groundwater. Because of these reasons, the plausibility check of the modeled nitrate concentrations in the leachate at the state level should not aim to assess the performance of the modeling at a specific location and time. Instead, consistency should be tested for larger regional reference units, such as groundwater bodies, with the goal of demonstrating that “hot spot” regions and “unproblematic” regions are equally well represented in the modeling.

Against this background, the good accordance of the modeled and observed nitrate concentrations prove, on the one hand, that the reactive N transport in the soil, i.e., the processes of N immobilization and denitrification, are adequately represented in terms of their magnitude and spatial distribution. On the other hand, it is thus demonstrated that the modeled nitrate concentration in the leachate represents a suitable reference value for assessing both the nitrogen reduction required to achieve the groundwater quality objective (50 mg/L) and the effects of concrete reduction measures on the development of nitrate concentrations in the groundwater.

The general question is whether the locations of the monitoring wells in the groundwater bodies are representative. This is because the nitrate concentrations observed in those locations can only be considered representative in this case. If, for example, a groundwater monitoring well was established in a cropland, but the cropland covers only a very small portion of the corresponding groundwater body, the representativeness of the monitoring well and the measured nitrate concentration, respectively, is not ensured. Therefore, the validation of the model results on the nitrate concentration in the leachate, especially at the state level, should be complemented in the future by an assessment of the positioning of the monitoring wells.

The denitrification in the groundwater simulated in the WEKU model was confirmed by the N₂/Ar measurements in the groundwater, so that both the hydrogeochemical milieu characterization for the regional designation of denitrification conditions in the groundwater and the assumed denitrification kinetics in the groundwater, as well as the modeled travel times in the groundwater, are adequately reproduced. In Hesse, however, the N₂/Ar measurements in the groundwater started just recently, i.e., there are no longer time series available. Accordingly, the conclusions so far are not statistically safe. The number of N₂/Ar monitoring sites and the frequency of the sampling will increase in the coming years, enhancing the reliability of the corresponding plausibility checks.

Due to a lack of monitoring data at the state scale to check the plausibility of the N inputs modeled separately for the individual input pathways, only the total sum of the modeled N inputs to the surface waters from diffuse and point sources could be validated. For this purpose, the measured N loads from 98 water quality monitoring stations were used, which showed a good agreement.

As with the modeled runoff, the good correspondence of the modeled total nitrogen inputs with the observed values at the surface water monitoring sites is no proof that the values determined for the individual input pathways are correct. However, there is no doubt that the contribution of the individual input pathways to the measured total N load is very different at the 98 gauge-based catchments used for the plausibility check (Figure 10). It is reasonable to assume that such a good agreement at so many gauges can only be achieved if the different contributions of the individual input pathways to the total load are well represented for each of the catchments. Accordingly, we conclude that the regional N emission pathways are sufficiently well represented in the model.

The manifold plausibility checks using the measured values in their combination give a good indication of the performance of the model system mGROWA-DENUZ-WEKU in Hesse. The confirmation of the model results based on the measured values leads to the overall conclusion that the N fluxes are adequately represented with this model system. Accordingly, it can be assumed that the effect of the changed N inputs, e.g., due to the reduced N fertilization in agriculture, can be predicted with some accuracy. The latter is especially important to confirm the usability of the model system for N management issues, e.g., for predictions on the effect of mitigation analyses.

8. Conclusions

With the model system mGROWA-DENUZ-WEKU, an instrument was created that considers the nitrate fluxes from their sources across all input pathways into groundwater and surface waters on a state-wide and spatially high-resolution basis. The intensive participation of the Hessian State Agency for Nature Conservation, Environment and Geology has ensured access to the best state-wide available databases, so that specific site characteristics are preserved in the model.

The multiple plausibility checks carried out showed that the results of the model system mGROWA-DENUZ-WEKU adequately represent the current nitrate pollution situation in the different Hessian landscapes. Accordingly, it can be expected that the effect of the changed N inputs can be predicted with a certain accuracy. The latter is especially important to underpin the applicability of the model system for N management issues, e.g., for predictions on the effect of mitigation analyses due to the reduced N fertilization in agriculture. For exactly this purpose, the model results are used in Hesse as a uniform and consistent information base for reporting and management plans in the EU WFD process.

In a more general sense, the plausibility checks showed that the model system is fundamentally suitable for representing the nitrate inputs to the water bodies on the state level. Although the availability of the model input data is limited in many countries, the potential for transferring the model system to other countries is possible. However, monitoring the data is required in any case for the calibration and validation of the models. As the study shows, it is only possible to carry out the plausibility checks in the manner described here if both the number of quantity and quality measuring points and the sampling frequency are sufficiently high, and if these data are made available centrally. Likewise, the verification of the model results for the nitrate concentration in the leachate can be improved by a higher measuring network density, as well as by the readiness of governmental institutions to provide the corresponding data.