Denitrification in Intrinsic and Specific Groundwater Vulnerability Assessment: A Review

Several groundwater vulnerability methodologies have been implemented throughout the years to face the increasing worldwide groundwater pollution, ranging from simple rating methodologies to complex numerical, statistical, and hybrid methods. Most of these methods have been used to evaluate groundwater vulnerability to nitrate, which is considered the major groundwater contaminant worldwide. Together with dilution, the degradation of nitrate via denitrification has been acknowledged as a process that can reduce reactive nitrogen mass loading rates in both deep and shallow aquifers. Thus, denitrification should be included in groundwater vulnerability studies and integrated into the various methodologies. This work reviewed the way in which denitrification has been considered within the vulnerability assessment methods and how it could increase the reliability of the overall results. Rating and statistical methods often disregard or indirectly incorporate denitrification, while numerical models make use of kinetic reactions that are able to quantify the spatial and temporal variations of denitrification rates. Nevertheless, the rating methods are still the most utilized, due to their linear structures, especially in watershed studies. More efforts should be paid in future studies to implement, calibrate, and validate user-friendly vulnerability assessment methods that are able to deal with denitrification capacity and rates at large spatial and


Introduction
Around the world, several water bodies do not meet the minimum quality standards necessary to ensure their utilization for anthropic activities [1][2][3] and groundwater quality is even declining over time [4]. This phenomenon is attributable to the presence of different pollutants in both surface and groundwater resources [5], due either to anthropogenic activities or natural processes [6]. Specifically, the subsurface pollution related to nitrogen surplus is considered a persistent and widespread issue, especially in agricultural areas [7], where the excess of nitrate (NO 3 − ) represents a major concern [8][9][10][11]. NO 3 can reduce the reactive nitrogen (Nr) levels in both deep and shallow aquifers [17,18]. The process occurs both in the saturated and unsaturated zones under mildly reducing conditions [19], with dissolved oxygen usually below 2 mg/L. It generally refers to the heterotrophic DNT, operated by facultative anaerobic bacteria, using a carbon source as electron donor and NO 3 − as electron acceptor to produce atmospheric nitrous oxide (N 2 O) and nitric oxide (NO), which can be further reduced to dinitrogen (N 2 ). Unfortunately, the DNT byproducts that are generated when DNT is incomplete are also considered greenhouse gases and ozone-depleting substances [20]. While heterotrophic DNT is the most widely known pathway to reduce NO 3 − in the subsurface, other processes can help to diminish NO 3 − concentrations, like dissimilatory NO 3 − reduction to ammonium [21]. This process is an anaerobic respiration of organic matter by chemoorganoheterotrophic microbes using NO 3 − as electron donor, but unlike DNT, it conserves bioavailable Nr in the system as ammonium. Other relevant processes ubiquitously occurring in subsurface environments are anammox [22] and autotrophic DNT via different inorganic electron donors [23], includign dissolved reduced species such as ferrous iron and sulfide or solid phases such as pyrrhotite and pyrite. Despite the many biogeochemical processes that can attenuate NO 3 − , fertilization often occurs at amounts exceeding both crop demand and the natural attenuation capacity of soils. This generates a NO 3 − surplus which can be readily transported by percolating waters through the soil profile and can promote NO 3 − accumulation in groundwater bodies, if oxic conditions prevail [24,25]. Moreover, a reliable identification of NO 3 − sources may be difficult due to the occurrence of multiple releasing sources and the complementarity of several geogenic and anthropogenic processes acting in the same area [26][27][28][29]. In addition, the occurrence of NO 3 − may suffer from consistent temporal variations [30], depending on: (i) the precipitation regime [31], (ii) the hydrogeological conditions [32], and (iii) the land use [33].
Several tools and techniques have been proposed to estimate groundwater pollution, ranging from analytical, graphical, and statistical methods. They have been greatly enhanced since the advent and the improvement of geographic information systems (GIS), which made the interpretation and representation of groundwater quality easier than before [34]. Among all available methodologies, groundwater vulnerability assessment tools have attracted researchers' attention due to their relatively straightforward applicability and the production of understandable results by non-specialists. The concept of aquifer vulnerability [35,36] aims to quantify the susceptibility of an aquifer to be adversely affected by any pollutant load imposed on the land surface. Groundwater vulnerability cannot be directly measured and quantified in the field since it represents a qualitative classification as a function of the main hydrogeological conditions of the study area, which could vary drastically from one place to another. The vulnerability can be classified in two types: (i) intrinsic and (ii) specific vulnerability [37]. The intrinsic vulnerability only depends on hydrogeological characteristics regardless the pollutant types. The specific vulnerability instead refers to the vulnerability of a given pollutant, involving in the evaluation processes both the pollutant's physical-chemical characteristics and its mass loading rates. In recent decades, different methodologies to assess groundwater vulnerability were proposed and tested around to world. They can be distinguished in four main categories: (i) deterministic or simulation models, (ii) overlay and index techniques, (iii) statistical tools, and (iv) hybrid methodologies [38]. All methodologies offer different advantages and, at the same time, they suffer from some drawbacks. For example, both statistical and deterministic models ask for a high-quality dataset and a massive quantity of information, especially for calibration and validation procedures, while overlay and indices methodologies could be affected by subjectivities in rates classification and weights assignment. Weighting and rating indices are easy-to-use methodologies, need relatively low amounts of data, are usually applicable to large domains, and can be easily modified and updated. On the other hand, they produce static snapshots of the aquifer vulnerability, which can be only used as a screening tool for regional and national scale management or planning. Conversely, process-based methods can produce a transient representation of groundwater quality, but they are mostly used on small domains (e.g., field or site scales) due to computational burdens and the requirement of high-quality and -frequency data.
In general, the groundwater pollution status is controlled by a wide range of physical and biogeochemical processes (e.g., dilution, sorption, and degradation), which contribute to decreasing pollutants concentrations. As a result, N loads of rivers and aquifer systems often reflect only a minor part of the N surplus within the corresponding river basin area due to decay processes during the aquifer transit and within the surface water bodies themselves [39]. Many studies have highlighted that among these decay processes, DNT has a large potential to reduce the relative vulnerability of water bodies [40,41]. Unfortunately, these attenuation processes are not always properly incorporated within the vulnerability assessment methodologies [42,43], but only considered indirectly in lumped parameters.
This paper deals with the last 20 years of literature on groundwater vulnerability and, specifically, it focuses on the introduction of NO 3 − attenuation processes within the evaluation assessment. The main difficulties in including, directly or indirectly, the DNT process in some methodologies, the improvements in results, and the effects on the consequent water management strategies are also discussed. With this review, we sought to highlight the main limitations and research gaps incurred in the past studies, delineating the future research needs aimed to obtain a reliable evaluation of groundwater intrinsic and specific vulnerability to NO 3 − , explicitly including DNT rates for both the saturated and unsaturated zones.

NO 3
− attenuation is controlled by several physical and biogeochemical processes occurring in the saturated and unsaturated zones. The understanding and quantification of these processes are mandatory to protect groundwater from the risk associated with point and non-point sources of NO 3 − pollution [44]. Figure 1 reports a schematic representation of the processes acting in the unsaturated and saturated zones.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 3 of 20 of groundwater quality, but they are mostly used on small domains (e.g., field or site scales) due to computational burdens and the requirement of high-quality and -frequency data.
In general, the groundwater pollution status is controlled by a wide range of physical and biogeochemical processes (e.g., dilution, sorption, and degradation), which contribute to decreasing pollutants concentrations. As a result, N loads of rivers and aquifer systems often reflect only a minor part of the N surplus within the corresponding river basin area due to decay processes during the aquifer transit and within the surface water bodies themselves [39]. Many studies have highlighted that among these decay processes, DNT has a large potential to reduce the relative vulnerability of water bodies [40,41]. Unfortunately, these attenuation processes are not always properly incorporated within the vulnerability assessment methodologies [42,43], but only considered indirectly in lumped parameters.
This paper deals with the last 20 years of literature on groundwater vulnerability and, specifically, it focuses on the introduction of NO3 − attenuation processes within the evaluation assessment. The main difficulties in including, directly or indirectly, the DNT process in some methodologies, the improvements in results, and the effects on the consequent water management strategies are also discussed. With this review, we sought to highlight the main limitations and research gaps incurred in the past studies, delineating the future research needs aimed to obtain a reliable evaluation of groundwater intrinsic and specific vulnerability to NO3 − , explicitly including DNT rates for both the saturated and unsaturated zones.

Denitrification in the Unsaturated and Saturated Zones
NO3 − attenuation is controlled by several physical and biogeochemical processes occurring in the saturated and unsaturated zones. The understanding and quantification of these processes are mandatory to protect groundwater from the risk associated with point and non-point sources of NO3 − pollution [44]. Figure 1 reports a schematic representation of the processes acting in the unsaturated and saturated zones. Among all biogeochemical processes, DNT is considered the predominant factor responsible for NO3 − attenuation in agricultural soil [45], representing an important pathway for Nr losses. Nevertheless, DNT also has a serious adverse environmental effect, Among all biogeochemical processes, DNT is considered the predominant factor responsible for NO 3 − attenuation in agricultural soil [45], representing an important pathway for Nr losses. Nevertheless, DNT also has a serious adverse environmental effect, being the principal source of N 2 O and NO, accounting for 70% of the N 2 O emitted annually from the biosphere into the atmosphere [46]. Several factors regulate DNT in the unsaturated and in the saturated zones; in fact, the DNT rate is a direct function of: (i) soil water content, (ii) hypoxic and anoxic conditions, (iii) NO 3 − concentrations, and (iv) organic carbon substrates, such as soil organic matter and crop residues [47]. Soil texture could also directly influence DNT, especially in the topsoil, via the formation of micro-niches where anoxic condition would prevail [48]. Nr mass loading rate significantly affects DNT, especially when the Nr supplied through fertilizer is higher than the crop Nr requirement [49,50]. Within the unsaturated profile, DNT could prevent leached NO 3 − from reaching the water table, if labile carbon sources are also present [51,52]; in any case, DNT rates could be significantly limited by the vertical unsaturated geochemical and microbiological variability.
DNT in the saturated zone has been less extensively studied with respect to DNT in the unsaturated zone. Despite this, DNT in groundwater systems is still relevant and mainly controlled by the same conditions as for the unsaturated zone. For instance, the presence of electron donors (organic carbon, reduced iron, and/or reduced sulfur) have been identified as the most limiting factor among with anoxic conditions, more than NO 3 − and nutrient availability [44]. The above-mentioned DNT controlling factors (Figure 1), both in the unsaturated and saturated zones, are usually considered in vulnerability assessment methodologies to indirectly account for DNT where direct measurements are not available. Furthermore, as will be further discussed in the following chapters, the main factor considered in most used methodologies (namely the rating and weights methods) is the soil texture, which is easier to measure and analyze compared to the other factors controlling DNT in the unsaturated or saturated zone.
Based on these considerations, Table 1 shows the DNT rates for the unsaturated zone, and more specifically in the topsoil, obtained from both laboratory and field studies, depending on the soil texture and for different land covers. The average DNT value for laboratory studies is 460 Kg-N/ha/year, with values that can reach up to 3600 Kg-N/ha/year, while for field studies the average DNT rate is 60 Kg-N/ha/year, with a maximum of 550 Kg-N/ha/year. It is evident that under controlled laboratory conditions, the DNT rates tend to be much higher (approximately one order of magnitude) than in field conditions, where local limiting factors may act to reduce the effectiveness of DNT both spatially and seasonally. Concerning the DNT rates in the saturated zone, it must be stressed that the studies available in the literature are not numerous, especially in the field, given the greater difficulty of the monitoring techniques required and the duration of the studies to be carried out. Nevertheless, recent studies [53,54] reported DNT rates up to 320-390 Kg-N/ha/year, underlining the fact that the contribution of DNT processes within the saturated zone cannot be neglected in the definition of the vulnerability of a given territory.

Materials and Methods
For this review, several studies focusing on intrinsic and specific groundwater vulnerability assessments performed in the last 20 years (2000-2020) worldwide were considered. The Scopus database (Elsevier) was probed as the main worldwide peer-reviewed manuscripts archive, using specific combinations of keywords, such as: groundwater vulnerability, DNT, nitrate vulnerable zone, nitrate attenuation factors, specific and intrinsic vulnerability. Given the large number of studies and applications focusing on the groundwater vulnerability topic, only the contributions published in peer-reviewed international journals were considered. Furthermore, the use of the keyword DNT in combination with at least others two keywords was necessary to exclude all the studies not related to the chosen topic. Specifically, all the studies on standard rating methods, such as DRASTIC, GOD, or SINTACS, without a clear reference to DNT were excluded, along with process-based models not adequately adapted to vulnerability and DNT assessment purposes. Similarly, this was done for the statistical and hybrid methodologies.
This review is not meant to serve as a general review on vulnerability assessment methods, since recent complete studies are already present in the literature [36,38,43], but aimed only to revise the role of DNT in groundwater vulnerability assessment studies; accordingly, only those papers where DNT was explicitly considered were retained. All the reviewed studies were divided into four main categories ( Figure 2) discussed separately: (i) rating methods, (ii) numerical models, (iii) statistical methods, and (iv) hybrid methodologies. This review offers an overall explanation of how already known vulnerability methodologies incorporate the process of DNT without discussing every single application, focusing the discussion on methodologies that showed an implicit or explicit involvement of DNT in the vulnerability assessment. Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 20

Results and Discussion
In this section, the main categories highlighted in Figure 2 will be discussed to assess the role of DNT in enhancing the already available intrinsic and specific groundwater vulnerability assessment methods.

Rating Methods for Nitrate Vulnerability Assessment
Among all methodologies to assess aquifer vulnerability, the rating methods are surely the most widely utilized [38]. The main advantage lies in their straightforward applicability and in their highly understandable results. The pioneer methods of this class are DRASTIC [77] and GOD [78], followed by many other methodologies specifically created to assess the intrinsic aquifer vulnerability, such as AVI [79], EPIK [80], SINTACS [81], SEEPAGE [82], and many more. A comprehensive review of these methodologies is available in Machiwal et al. [38]. These methods have been widely applied both on watershed, regional, and even global scales, producing reliable results [83][84][85][86][87]. They consist of the combination of thematic layers, representing the different physiographic attributes that control the groundwater vulnerability to pollution of a given area; the most utilized factors are: (i) depth to the water table, (ii) recharge, (iii) unsaturated and saturated media, (iv) soil texture, (v) topographic slope, (vi) hydraulic conductivity, (vii) hydraulic resistance, and (viii) presence of phreatic or confined aquifers. Each parameter is classified using a range of rates (generally 1 to 10) and multiplied for a specific weight, indicating the individual parameter importance and magnitude. These lumped parameters incorporate most physical attenuation processes (such as dilution and travel time) but disregard biogeochemical processes (such as DNT) that may occur along the pathway between the pollution sources (usually at land surface) and the considered target, namely the saturated portion of the aquifer and/or the surface water bodies towards which groundwater discharges ( Figure 3).

Results and Discussion
In this section, the main categories highlighted in Figure 2 will be discussed to assess the role of DNT in enhancing the already available intrinsic and specific groundwater vulnerability assessment methods.

Rating Methods for Nitrate Vulnerability Assessment
Among all methodologies to assess aquifer vulnerability, the rating methods are surely the most widely utilized [38]. The main advantage lies in their straightforward applicability and in their highly understandable results. The pioneer methods of this class are DRASTIC [77] and GOD [78], followed by many other methodologies specifically created to assess the intrinsic aquifer vulnerability, such as AVI [79], EPIK [80], SINTACS [81], SEEPAGE [82], and many more. A comprehensive review of these methodologies is available in Machiwal et al. [38]. These methods have been widely applied both on watershed, regional, and even global scales, producing reliable results [83][84][85][86][87]. They consist of the combination of thematic layers, representing the different physiographic attributes that control the groundwater vulnerability to pollution of a given area; the most utilized factors are: (i) depth to the water table, (ii) recharge, (iii) unsaturated and saturated media, (iv) soil texture, (v) topographic slope, (vi) hydraulic conductivity, (vii) hydraulic resistance, and (viii) presence of phreatic or confined aquifers. Each parameter is classified using a range of rates (generally 1 to 10) and multiplied for a specific weight, indicating the individual parameter importance and magnitude. These lumped parameters incorporate most physical attenuation processes (such as dilution and travel time) but disregard biogeochemical processes (such as DNT) that may occur along the pathway between the pollution sources (usually at land surface) and the considered target, namely the saturated portion of the aquifer and/or the surface water bodies towards which groundwater discharges (Figure 3). It should be mentioned that some biogeochemical processes are often a function of soil properties (such as texture and water content), therefore, many end users of rating methods claim that natural attenuation processes are accounted for, even though not ex- It should be mentioned that some biogeochemical processes are often a function of soil properties (such as texture and water content), therefore, many end users of rating methods claim that natural attenuation processes are accounted for, even though not explicitly. Specifically, soil texture could influence the redox values suitable for DNT, which were found to be more favorable in clayey textured soils than in sandy textured soils [88,89]. Nevertheless, other important factors influencing DNT, like organic carbon substrates and Nr mass loading rates, are generally completely neglected. The same assumption is made for the unsaturated and saturated zone where the pollutant attenuation can be regulated by several processes that depend on their physical and hydraulic properties (e.g., thickness, hydraulic conductivity, and presence of organic matter). Both saturated and unsaturated zones are commonly classified considering only the lithological characteristic and the thickness of the media. Moreover, all methods offer a large range of rates, especially for saturated and unsaturated media classification, implying an inherent subjectivity (expert judgment). In any case, the classification is usually made considering dispersion and diffusion capacity, neglecting the main biogeochemical factors that can occur within the saturated and unsaturated zones. The SEEPAGE methodology [82] tried to overcome this problem, introducing a potential attenuation parameter referring to the soil and subsoil characteristics. Specifically, the introduced attenuation parameters based its classification on several factors that can be easily related to DNT, such as: (i) soil surface and subsoil texture, (ii) topsoil pH, (iii) organic matter content of topsoil, (iv) soil drainage class, and (v) soil permeability. Unfortunately, the methodology has not been so widely applied in comparison to DRASTIC and SINTACS methods and only a few applications have been documented so far [90,91].
On the other hand, much more attention has been paid to the determination of the different NO 3 − sources, so much so that there are numerous studies that have correlated the rating methods to land use [92][93][94] and to probabilistic models [95,96].

Nitrate Vulnerable Zones
A step forward has been made with the introduction of several methodologies aimed at evaluating nitrate vulnerable zones (NVZ). Conserving the same structure of a rating method, specific parameters connected to NO 3 − leaching and attenuation were introduced in the evaluation procedure. An example is the agricultural NO 3 − hazard index named IPNOA, which was proposed by Padovani and Trevisan [97]. IPNOA facilitates the assessment of hazard factors contributing to NO 3 − load and control factors able to amplify or decrease the NO 3 − load. Specifically, the hazard factors are chosen among agricultural sources of Nr input which could have a potential impact on groundwater resources, such as: (i) applied organic fertilizers, (ii) synthetic fertilizers, and (iii) wastewaters, while the control factors used to evaluate NO 3 − behavior in site-specific conditions and farming practices are: (i) Nr content in soil, (ii) climate characteristic, (iii) agronomic practices, and (iv) irrigation techniques. For example, tillage operations could greatly influence DNT [98,99], as well as the soil Nr content and the fertilizer application rate [47]. The methodology was applied by several authors [33,[100][101][102] and was also coupled with classical rating methods, such as SINTACS, with the aim of improving NO 3 − vulnerability assessment [103,104]. Another approach for NVZ delineation is the calculation of N budgets and surpluses on the farmland scale [105], watershed scale [106][107][108][109][110], or even at the subcontinental scale [111][112][113]. It consists of differentiating between the soil Nr inputs and outputs, where inputs usually are: (i) synthetic fertilizers, (ii) sewage sludge and animal manure, (iii) biological N fixation, and iv) atmospheric deposition, while the outputs are: (i) the amount of crop harvested and removed from the field, (ii) crop residues, (iii) leached Nr, and (iv) degassed N. This methodology offers a mappable output of the available Nr content for leaching. Rebolledo et al. [114] presented an overlay index method for NVZ identification through the analysis of three groups of factors: the first was related to intrinsic vulnerability (water table, saturated hydraulic conductivity, and aquifer type), the second considered the attenuation factors acting in the soil and specifically related to DNT capacity (organic carbon, pH, and texture), and the third indicated the pollutant infiltration capacity (slope, precipitation, and evapotranspiration). Rates and weights were then assigned using Pearson's r correlation and analytical hierarchical process (AHP) to limit subjectivity. Arauzo et al. [115] combined the N-surplus approach with a rating methodology specially developed for NVZ. Similarly, Cameira et al. [116] developed a global risk index (GRI), integrating the N budget with water surplus (irrigation and rain), soil permeability, unsaturated zone residence time, and surface morphology. In these methodologies, DNT is related to different parameters such as soil texture, soil pH, Nr input, and unsaturated zone residence time. Indeed, DNT is more likely to occur when the travel time through the unsaturated zone is long [44]. Following the same concept, Busico et al. [117] proposed a rating methodology for NVZ using all those parameters directly connected to NO 3 − concentration in groundwater, combining the classical parameters of the unsaturated zone media, surface morphology, and recharge with the natural soil Nr content and the amount of fertilizers recommended for different kind of crops. The methodologies described so far have been applied in several case studies and compared to canonical vulnerability assessment methods and were found to deliver more realistic results [117][118][119].

Statistical Methods
The statistical methodologies aim to assess groundwater vulnerability to NO 3 − mapping pollution probabilities at a specific site. The computation is based on the correlations between some parameters describing the aquifer properties and the sources and occurrences of pollutants, such as: (i) aquifer type (lithology or confinement condition), (ii) Nr mass loading rates, (iii) water physical-chemical variables (major ions, trace elements, and isotopic composition), and (iv) soil characteristics (texture, organic carbon, and pH). To fulfill this purpose, these methodologies require extensive datasets on both physiographic and geochemical properties of the study area. They may incorporate a wide range of inputs related to both the anthropogenic stresses and the natural conditions of the groundwater system [120]. Among the various statistical tools, the most utilized for aquifer vulnerability assessment to NO 3 − are: (i) logistic regression [121], (ii) multiple linear regression [122], and (iii) geostatistical techniques [123]. Logistic regression allows for the evaluation of those natural and anthropogenic factors that could impact the probability that NO 3 − concentrations will exceed a specified threshold. Commonly, statistical methods involve numerous factors related to the attenuation processes and to DNT, like soil characteristics and Nr loads [124][125][126][127][128][129]. Moreover, some authors identified the anoxic condition as the main attenuation factor [126,[130][131][132]. Boy-Roura et al. [133] built a multiple linear regression model using five main variables, one of which was identified as indicator of DNT occurrence. The variable was assessed using N-stable isotopes in NO 3 − and dissolved oxygen.

Numerical Models
With respect to the qualitative methods (rating methods), the numerically and physically based methods for the assessment of aquifer vulnerability to NO 3 − involve the analysis of all the relevant processes occurring within the unsaturated and the saturated systems. Numerical models subdivide the tridimensional domain into cells or elements (depending on the method used) where the flow, dispersive transport, and reactions among each cell are computed by differential equations of flux and mass conservation solved iteratively, starting from user-defined initial and boundary conditions. Thus, with numerical methods, an accurate temporal and spatial evolution of the fate and transport of NO 3 − in the subsurface can be achieved. In fact, these methodologies aim to determine the specific vulnerability for diffuse pollution sources rather than the intrinsic vulnerability of an aquifer, emphasizing the physical, chemical, and biological processes controlling the fate and transport of contaminants in both the unsaturated and saturated zones. Moreover, process-based models also provide the temporal variation of vulnerabil-ity, a concept often neglected by other methods. Among all the numerical models, the most used are: GLEAMS [134], SUTRA [135], MODFLOW [136] coupled with MT3DMS [137], FEFLOW [138], PHT3D [139], HYDRUS [140,141], and SWAT [142]. All of these models explicitly involve Nr attenuation processes, such as DNT, in their modules. For example, the N module of SWAT intrinsically considers various NO 3 − reduction processes, such as mineralization, volatilization, and DNT, which is simulated in the soil layer as a function of the amount of NO 3 − , organic carbon, and temperature. The sink and source package of MT3DMS was designed to simulate chemical processes, such as NO 3 − decay, following a first-order degradation reaction or fully integrating the reaction stoichiometries of the main terminal electron acceptors processes [143]. GLEAMS incorporates a N cycling algorithm, which considers N plant uptake and fixation among with fertilizers and manure applications. Due to their high versatility, despite these models are recognized as valuable tools for groundwater vulnerability assessment, they find application in a large variety of settings, such as: pollution transport in groundwater and surface waters, groundwater balance and depletion, seawater intrusion, and more.
This picture shows that groundwater vulnerability assessment based on contaminant transport modeling can provide more quantitative evaluation than previously discussed rating methods. On the other hand, the model's large number of input parameters and the uncertainty associated with them may be considerable issues. In fact, to implement such simulations, an accurate characterization of all flow and transport processes acting at a given site is a prerequisite, and, in addition, there is also the need to deeply understand the geochemical environment [144].
Despite this, the possible applications of numerical models are quite numerous. Some examples of aquifer vulnerability to NO 3 − measured using MODFLOW/MT3DMS models are found in USA [145], Europe [146], and Africa [147]. Uhan et al. [148], assessed NO 3 − groundwater vulnerability for an alluvial aquifer combining the output of three different numerical models: (i) groundwater recharge (GROWA), (ii) NO 3 − leached from the soil profile (SWAT), and (iii) groundwater flow velocities (FEFLOW). Cui et al. [149] assessed the whole Florida shallow aquifer vulnerability through the Nr removal and transport rate in the unsaturated zone. Hansen et al. [150] developed the site-specific concept for aquifer nitrate vulnerability assessment (SCANVA) through the combination of a 3D geological reconstruction, a groundwater simulation based on MODFLOW-2000, and a complete hydro-geochemical assessment. Huan et al. [151][152][153] proposed a specific vulnerability assessment using HYDRUS-1D model and FEFLOW, emphasizing the role of DNT within the unsaturated zone using an enhanced convection-dispersion equation. The Nr transfer time and NO 3 − susceptibility in pumping wells were finally used to visualize the specific vulnerability.

Hybrid Methodologies
Some authors tried to enhance the results of the assessment of vulnerability to NO 3 − , combining classical rating methods along with process-based and statistical elaborations. Sometimes, these new hybridizations introduced new parameters related to DNT that are usually not considered in the standard methodologies. Keuskamp et al. [68] hybridized N budget with a basin-scale process-based model to simulate the fate of Nr in soils and groundwater, along with N 2 O production at the European scale. The model offers the output calculation of travel time, which integrates in its computation the DNT rate for the main European countries. Aschonitis et al. [154,155] developed a simplistic N-budget approach described by indices named LOS, which account for DNT, NH 3 volatilization, mineralization, and nitrification. The new indices involve detailed information on landscape management and climate and are calibrated using the results of a GLEAMS model thought regression analysis. Kazakis and Voudouris [156], following the example of Huan et al. [157], who subjectivized the DRASTIC framework, replaced some parameters in the DRASTIC-PAN methodology. This new method, integrated within the LOS indices, indicates the soil Nr losses based on climatic, soil, and topographic data. Similarly, Bu-sico et al. [11,158] hybridized the SINTACS method, called SINTACS-SVN and SVAP, introducing the LOS indices inside the evaluation process to assess intrinsic and specific groundwater vulnerability to NO 3 − and SO 4 2− pollution. In both methodologies, the weights and scores were calibrated using observed NO 3 − concentrations and pollution factors following the routine of canonical process-based methodologies, making the evaluation less subjective than standard rating methods. These indices proved to produce reliable results in different hydrogeological settings [159]. Jia et al. [160] developed the DRANTHVP method, optimizing weights and rates of the DRASTIC method by using the statistical analysis of projection pursuit dynamic clustering (PPDC). Moreover, during the choice of parameters, the authors introduced the NO 3 − attenuation factor to replace the soil media factor. The proposed classification is based on the difference of soil moisture content, porosity, and organic matter content, factors that are directly linked to DNT rate [161,162]. Vouthckova et al. [163] integrated subsurface redox conditions inside the DRASTIC classification. This parameter is considered indicative of DNT as the NO 3 − removal by microbial oxidation depends on the availability of reductants species, such as pyrite and organic matter [164].

Concluding Remarks
The assessment of groundwater vulnerability to NO 3 − is a valuable tool for the implementation of sustainable groundwater management plans. However, the standard intrinsic vulnerability assessment methods can lead to unreliable results, since they do not consider all NO 3 − attenuation processes occurring in both the unsaturated and saturated zones, as proved by validation procedures using observed NO 3 − concentrations providing feedback on the reliability of the applied methodology. Accordingly, in the last 20 years, the classical rating methodologies have been deeply modified to account for the pollutant dilution and attenuation processes, with the aim of enhancing their reliability. Still, the analysis of all the studies considered in this review has shown that almost 50% of groundwater vulnerability assessments are performed using rating methodologies, followed by statistical methods, numerical models, and hybrid methods (Figure 4).  On the other hand, in the statistical and numerical models, the concept of validation is intrinsically linked to the model structures. Numerical models require validation on multiple steps and are always preceded by calibration procedures aimed at quantifying the model's ability to simulate the real system. The statistical methods follow the same concept of calibration and validation procedures, where the spatial distribution of the socalled explanatory variable (aquifer type, soil characteristics, agricultural Nr load, and land cover) is discretized according to water quality. Moreover, statistical methods always consider observed data (pollutant concentration, dissolved oxygen, soil pH, and more) in the evaluation process. Unfortunately, due to their complex structure and necessity of The reason for this result is likely due to the low data requirement and simplicity of elaboration needed for the rating methodologies with respect to the other methods. Inside the group of rating methodologies (Figure 3), the NVZ applications introduced some parameters directly connected to DNT, such as soil organic carbon, agricultural practices, and Nr input, but still retained the structure of a classical rating method. Indeed, in the NVZ methods, most of the studies used a spatial representation of the validation procedure, projecting the NO 3 − concentrations over the vulnerability classification, but without accounting for groundwater travel times.
On the other hand, in the statistical and numerical models, the concept of validation is intrinsically linked to the model structures. Numerical models require validation on multiple steps and are always preceded by calibration procedures aimed at quantifying the model's ability to simulate the real system. The statistical methods follow the same concept of calibration and validation procedures, where the spatial distribution of the so-called explanatory variable (aquifer type, soil characteristics, agricultural Nr load, and land cover) is discretized according to water quality. Moreover, statistical methods always consider observed data (pollutant concentration, dissolved oxygen, soil pH, and more) in the evaluation process. Unfortunately, due to their complex structure and necessity of specific data, they are mainly confined to site-specific assessments, even if some examples of regional and continental applications exist. Similarly, the hybrid methods suffer from a complex structure, often based on regression analysis and statistical correlation, demanding more data and sometimes requiring specific outputs produced by complex numerical models. Despite the increased complexity, the hybrid methodologies showed a higher reliability, confirmed by the validation procedure performed in all the studies.

The Role of Denitrification in Groundwater Vulnerability Assessment
In summary, it is possible to state that the involvement of NO 3 − decay processes in the evaluation of groundwater vulnerability assessment, and specifically of DNT, directly increased the final reliability of the vulnerability evaluation. This is particularly evident in rating applications, where the methodologies that lumped the parameters directly connected with DNT showed better results than standard DRASTIC, SINTACS, AVI, or GOD methods [11,117,156,158,160]. However, this improvement pays the price for an increase in model complexity and data requirement, a common drawback also for numerical and statistical models.
Thanks to the recent increase in the availability of online and free datasets, required by open science and open data protocols, such drawbacks could be overcome in future studies. For example, Keuskamp et al. [68] provided an average yearly DNT rate on the river basin scale, which can be easily incorporated in statistical and numerical models. Also, Kumar et al. [7] proposed a NO 3 − leaching assessment on the European scale. The results of these computations, along with many others, could be eventually used to enhance site assessment vulnerability to NO 3 − , even in those areas where a direct measurement of DNT rates is not available. Furthermore, where data are scant, new approaches to evaluate groundwater vulnerability have been recently proposed [165] that leave behind the concept of rating methods to consider NO 3 − contamination, suggesting specific sub-indicators of exposure, sensitivity, and adaptive capacity according to the indications of the Intergovernmental Panel on Climate Change [166,167].
Summarizing, standard vulnerability methods, not considering all the biogeochemical processes occurring both in the saturated and unsaturated zones, still represent precious screening tools for a qualitative assessment of aquifer vulnerability especially at large scales, but when there is a need to gain quantitative, reliable, and time-dependent vulnerability assessments, the employment of more complex data-and time-requiring methodologies must be preferred.

Future Trends in the Application of Groundwater Vulnerability Assessment
The economic and social significance of a good land and water management is an overarching goal in the UN 2030 Agenda for sustainable development. Fundamental changes in the way societies produce and use natural resources are indispensable for achieving the sustainable development goals (SDGs) within the carrying capacity of ecosystems that provide essential services, which are critical for human well-being.
For this to happen, good governance and safeguarding of land and water resources should be implemented in an integrated manner, involving all stakeholders and building monitoring and assessment tools based on sound, scientifically formulated knowledge. For instance, boosting yields and protecting crops to assure food production has often been achieved via the extensive use of agrochemicals (fertilizers, pesticides, etc.), with the consequence that many aquifers in heavily anthropized lowlands record high values of these contaminants and, thus, have recently been classified as NVZ. To adequately address these issues, decision makers must look at the latest scientific knowledge, which can support productive food systems (e.g., SDGs 2 and 12) through sustainable soil, land, water, nutrient, and pest management. Refining the groundwater vulnerability assessment methods, including DNT in the evaluation process, goes precisely in this direction.
In this context, future studies should be carried out focusing on the 'source-pathwaytarget' approach [168], where source indicates potential contamination due to land-use activities, pathway refers to the fate and transport of contaminants from the source to the target, and the target is the environment receiving the contamination load. In the specific case of groundwater vulnerability, the pathway should be identified in the movement of the infiltrating water across the unsaturated zone and within the saturated zone, while the target could be identified in the aquifer itself and in the groundwater-dependent ecosystems (GDEs), which could be endangered by the high concentration of pollutants discharging from the aquifer to the surface water bodies. It is known that the concentration of pollutants depends on the characteristics of the source (loading rate and duration) but also on the characteristics of the pathway where dilution, dispersion, and/or attenuation may occur. For this specific reason, it is fundamental to include an explicit evaluation of DNT rates for the saturated and unsaturated zones to account for the attenuation capacity of the pathway [11,42]. Although the monitoring and characterization of the unsaturated zone is complex and time consuming, future studies need to fulfill this task to gain more realistic estimates in groundwater vulnerability assessments [38]. Moreover, since the processes occurring along the pathway are also highly variable in time and space [30], there is a need for improved understanding of DNT rate variability over the hydrologic year, both as a function of the natural seasonality of the attenuation processes and as a function of the different agricultural practices that are carried out at specific times of the year (fertilization, plowing, harvesting, etc.), which, in turn, may affect the factors controlling DNT in the unsaturated/saturated zones. Based on this information, it will be necessary to implement the groundwater vulnerability assessment tools with approaches and techniques that are able to account for this seasonality (like numerical models). This improvement will lead to great benefits, especially in the management of land and water resources in GDEs [169]. In fact, GDEs are often located in flat coastal areas (e.g., wetlands, swamps, riverine plains, etc.) that are also ideal environments for agricultural activities [170]. This combination makes most of these habitats prone to adverse effects of water pollution induced by human activities. One of the most common examples is the establishment of eutrophication conditions due to excessive N load that may lead to the death of aquatic species. For this reason, the systematic mapping of groundwater vulnerability, accounting also for the fate of N species (thus including DNT), must be promoted to develop land uses and human activities according to the groundwater vulnerability, responding to important SDGs of the 2030 Agenda (e.g., SDGs 14 and 15).
Finally, the recent tendency to improve the reliability of the methods considered in this review and the definition of new approaches for groundwater vulnerability assessment demonstrates that a sound tool to properly manage water resources is considered necessary and urgent, not only by the wider research community, but even by decision makers and land and water managers. In particular, the most coveted result is to provide a userfriendly tool based on the latest scientific knowledge, which would allow the definition of management plans capable of considering the most effective allocation of mitigation strategies and achieving the highest reduction of adverse impacts within a relatively short time and available budgets.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the Scopus database restriction access only to subscribers.

Conflicts of Interest:
The authors declare no conflict of interest.