Among all methodologies to assess aquifer vulnerability, the rating methods are the most widely utilized [38]. The main advantage lies in their straightforward applicability and in their highly understandable results. The pioneer methods are DRASTIC [77] and GOD [78], followed by other methodologies such as AVI [79], EPIK [80], SINTACS [81], and SEEPAGE [82]. A comprehensive review of these methodologies is available in Machiwal et al. [38]. These methods have been widely applied both on watershed, regional, and global scales [83,84,85,86,87]. They consist of the combination of thematic layers, representing the physiographic attributes that control groundwater vulnerability to pollution; 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 and multiplied for a specific weight, indicating the individual parameter importance and magnitude. These lumped parameters incorporate most physical attenuation processes, such as dilution, but disregard biogeochemical processes, such as denitrification (DNT) (Figure 1).

Some biogeochemical processes are often a function of soil properties, 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 are more favorable in clayey textured soils than in sandy textured soils [88,89]. Nevertheless, other important factors influencing DNT, like organic carbon substrates and nitrogen mass loading rates, are generally completely neglected. The SEEPAGE methodology [82] tried to overcome this problem, introducing a potential attenuation parameter referring to the soil and subsoil characteristics [90,91]. A step forward has been made with the introduction of methodologies aimed at evaluating nitrate vulnerable zones (NVZ). Conserving the same structure of a rating method, specific parameters connected to NO₃⁻ leaching and attenuation were introduced in the evaluation procedure. An example is the agricultural NO₃⁻ hazard index IPNOA [97] which defines hazard factors contributing to NO₃⁻ load and control factors able to amplify or decrease the NO₃⁻ load. 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. These methodologies 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].

The statistical methodologies aim to assess groundwater vulnerability to NO₃⁻ 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). Among the statistical tools, the most utilized for aquifer vulnerability assessment to NO₃⁻ are: (i) logistic regression ^[1], (ii) multiple linear regression ^[2], and (iii) geostatistical techniques ^[3]. Commonly, statistical methods involve numerous factors related to the attenuation processes and to DNT, like soil characteristics and Nr loads ^{[4][5][6][7][8][9]}. Moreover, some authors identified the anoxic condition as the main attenuation factor ^{[6][10][11][12]}. 

With respect to the qualitative methods (rating methods), the numerically and physically based methods for the assessment of aquifer vulnerability to NO₃⁻ 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 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, numerical methods aim to determine the specific vulnerability for diffuse pollution sources rather than the intrinsic vulnerability of an aquifer, and to provide the temporal variation of vulnerability, a concept often neglected by other methods. Among numerical models, the most used are: GLEAMS ^[13], SUTRA ^[14], MODFLOW ^[15] coupled with MT3DMS ^[16], FEFLOW ^[17], PHT3D ^[18], HYDRUS ^[19][20], and SWAT ^[21]. All of these models explicitly involve Nr attenuation processes, such as DNT, in their modules and 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 a considerable issues. Despite this, the possible applications of numerical models ^[22][23][24] and their combinations ^[25][26] are quite numerous . 

Some authors tried to enhance the results of the assessment of vulnerability to NO₃⁻, combining classical rating methods with process-based and statistical elaborations. Sometimes, these hybridizations introduced new parameters related to DNT that are usually not considered in the standard methodologies. Keuskamp et al. ^[27] hybridized N budget with a basin-scale process-based model to simulate the fate of Nr in soils and groundwater. Aschonitis et al. ^[28][29] developed a simplistic N-budget approach described by indices named LOS, which account for DNT, NH₃ volatilization, mineralization, and nitrification. Busico et al. ^[30][31] hybridized the SINTACS method, called SINTACS-SVN and SVAP, introducing the LOS indices in the evaluation process to assess intrinsic and specific groundwater vulnerability to NO₃⁻ and SO₄²⁻ pollution. Jia et al. ^[32] developed the DRANTHVP method, optimizing weights and rates of the DRASTIC method by using the statistical analysis of projection pursuit dynamic clustering (PPDC). Vouthckova et al. ^[33] integrated subsurface redox conditions in the DRASTIC classification to account for DNT, as the NO₃⁻ removal by microbial oxidation depends on the availability of reductants species, such as pyrite and organic matter ^[34].