Development of a groundwater contamination index based on the agricultural hazard and aquifer vulnerability: Application to Portugal

Bycombiningthehazardwithamulti-parametergroundwa-tervulnerability,aspatiallyexplicitgroundwatercontaminationrisk,developedformainlandPortugal,wascom-putedfor1999and2009.Resultsshowanincreasefrom8,800to82,679haoftheterritoryratedwithaveryhigh contamination risk. The priority areas were successfully screened by the Index, coinciding with the current Vul-nerableZones,althoughadditional hotspots were detected insouthernPortugal.Percolation, including both irrigation activity and precipitation, was found to be a key driver for the groundwater contamination risk due to its opposite effects in the hazard and in the vulnerability. Reducing nitrogen leaching may be insuf ﬁ cient to reduce the risk of nitrate contamination ifthere is a relatively larger reduction in precipitation. This indexis particularly useful when appliedtocontrasting situationsofvulnerability and hazard,whichrequiredistinctmitigation measures to mitigate groundwater contamination. © 2021


Introduction
The Mediterranean climatic regions have been identified as a hotspot for climate change (de Sherbinin, 2014;Diffenbaugh and Giorgi, 2012;Giorgi and Lionello, 2008).Prolonged dry periods, higher average temperatures and decrease in precipitation (Founda et al., 2019;Linares et al., 2020) may cause a 30-50% decline in freshwater availability (Milano et al., 2013).Coupled with increased anthropogenic pressure and higher water demand predicted for this region (Cramer et al., 2018), the existing groundwater overexploitation can be further aggravated (e.g.Goebel et al., 2019) and increase groundwater vulnerability (Nistor, 2019).
Agriculture is considered one of the major sources of diffuse groundwater pollution in Europe (Harrison et al., 2019) through the overapplication of agricultural nitrogen (N) inputs.This can lead to nitrate leaching to surface-and groundwater (Wang et al., 2019).High nitrate concentrations in drinking water are associated with several public health issues like colorectal (Schullehner et al., 2018), bladder, kidney and brain cancer (Ward et al., 2018) and methemoglobinemia (Sanchez-Echaniz et al., 2001).The EU Nitrates Directive 91/676/EEC (ND) stipulates a maximum value for human consumption of 50 mg NO 3 − L −1 .The imposition of a legal threshold in nitrate concentration in drinking water differentiates between nitrate contamination and pollution, depending on whether it is below or higher than the threshold, respectively.
The gross nitrogen balance (GNB) is an useful agri-environmental indicator of potential environmental N losses to the air, soil and water (Leip et al., 2011).Despite its simplicity, several authors have validated its usefulness in predicting nitrate leaching (Blicher-Mathiesen et al., 2014;Dalgaard et al., 2012;Wick et al., 2012).Furthermore, sitespecific leaching fractions have been applied to the GNB in order to estimate nitrate loads from below the root zone to deeper groundwater (e.g., De Vries et al., 2011).While this approach provides a proxy for the groundwater contamination hazard, it often fails to include the intrinsic characteristics of groundwater that influence the transport to the phreatic zone (i.e., vulnerability).
Another approach is to apply groundwater vulnerability indexes (e.g., DRASTIC; USEPA, 1985) to years where detailed data are available.However, these models focus on groundwater vulnerability assessment rather than on the contamination hazard and groundwater contamination risk.Thus, they can indicate high vulnerability, but no pollution risk given the absence of a contamination load (Cameira et al., 2021).By integrating the hazard and vulnerability concepts, the present paper aims to present a more holistic approach to assess and to rate the groundwater N contamination risk at a national scale.This methodology allows the identification of the main risk factor in each aquifer and therefore to assist technicians and policymakers in targeting the mitigation measures.To accomplish that, a risk-based index was determined for Portugal, using extensive data for the years 1999 and 2009, by coupling (i) an hazard module integrating the annual water balance surplus (precipitation + irrigation) and the N leaching obtained from the disaggregation of the previously calculated N gross balance into the different environmental losses, and (ii) a vulnerability module, that accounts for soil texture, land slope and nitrate residence time in the vadose zone.

Materials and methods
The methodology was applied on a national scale, considering mainland Portugal.The territory presents a considerable heterogeneity in climate, geology and land use in a relatively small area (92,145 km 2 ).

Precipitation in mainland Portugal: spatial and temporal patterns
Portugal is located in the transitional region between the sub-tropical anticyclone and sub-polar depression, with latitudes and longitudes ranging from 36°56′ and 42°09′N and 6°10 and 9°34′W, respectively.It is characterized as having mild and rainy winters, and warm and dry summers, the typical Mediterranean climate, though there are considerable climatic spatiotemporal variability due to different orographic conditions, particularly precipitation (Santos et al., 2017).The average annual precipitation for the historical period of 1970-2000 was 882 ± 172 mm yr −1 (Belo-Pereira et al., 2011).Whereas the mountainous regions in Northern Portugal are relatively rainy throughout the year (average = 2000 mm yr −1 ; e.g., NW Portugal), most of the precipitation occurs during winter in the flatter Southern Portugal (mean ≈500 mm yr −1 ).Portela et al. (2020) have analyzed the long-term rainfall trend in the latter regions  and found a marked decline.Conversely, Portela et al. (2014) showed that rainfall changes in mainland Portugal did not follow any significant trend, except from a considerable decline in March.Nonetheless, according to Belo-Pereira et al. (2011), precipitation decreased about 9 and 19% in 1999and 2009compared to the average historical data (1960-2000), particularly in NW Portugal (45%; Fig. 1).

Estimating leaching from nitrogen balances
The GNB previously calculated by Serra et al. (2019) at the municipality scale for the years 1999 and 2009 was used to estimate nitrate (NO 3 − ) leaching.The authors took a mass balance approach between agricultural N inputs and outputs (Eq.( 1)), where the former includes Fig. 1.Precipitation for 1999 and 2009 (SNIRH, 2019) and the spatial distribution of the main aquifer systems and the four hydrogeological units as well as the hydraulic conductivity (Tóth et al., 2014).
atmospheric N deposition (N dep ), biological N fixation (N bnf ), gross manure N (N man ), sludge (N ss ) and inorganic fertilisers applications (N fert ).The outputs include the N in crop (N crop ) and fodder production (N fodder ), crop residues removed from the field (N res ) or burnt in situ (N burnt ).All parameters are given in kg N ha −1 yr −1 .
Assuming no changes in the soil N stock, the GNB encompasses losses to the atmospheric and aquatic environments in the form of gaseous emissions of ammonia (NH 3 ), nitrous oxide (N 2 O), nitrogen oxides (NO x ) and molecular N (N 2 ), and N losses through surface runoff and leaching to ground-and surface water.By accounting these losses, the GNB can be disaggregated into different N balance typologies (Leip et al., 2011), which vary in purpose and data requirements.Here, we calculated the soil surface nitrogen balance (SSNB; kg N ha −1 yr −1 ) since Van Grinsven et al. (2012) found this a particularly useful proxy for N losses to ground-and surface water (e.g., dams, rivers) (Eq.( 2)).
where N gas is the sum of NH 3 , N 2 O and NO x emissions from manure management systems and agricultural soils (kg N ha −1 yr −1 ), N runoff is the N in surface runoff losses (kg N ha −1 yr −1 ), N leach is the nitrate leaching from below the root zone (kg N ha −1 yr −1 ) and N denit is the denitrification to molecular nitrogen (kg N ha −1 yr −1 ).If negative N balances occur, N leach and N denit were set to 0. These situations were attributed to statistical biases that affected the distribution of inorganic fertilisers in 1999 (Serra et al., 2019) and/or soil N depletion (Moklyachuk et al., 2019;Özbek and Leip, 2015).Table 1 summarizes data inputs and sources as well as the methods employed for each environmental N loss.Gaseous N losses as NH 3 , N 2 O and NO x were computed according to the manure management (animal housing and manure storage) and field application of fertilisers (including grazing).The IPCC (2019) emission factor (EF; %N-inputs) of 1% was used for direct N 2 O emissions, but for the Mediterranean regions of Portugal, the calculations were refined by using Cayuela et al. (2017) EFs for different irrigation systems (rainfed, sprinkler, surface and drip).Data on irrigated areas per crop and per irrigation systems at the municipality scale were available for 2009, while for 1999 the existing data included only irrigated areas at the agrarian region scale (Portugal Statistics, 1999, 2009).Consequently, data for 1999 were downscaled from the agrarian region to the municipality scale using Eq.(3) adapted from Cameira et al. (2019).
where IA m99 and IA m09 are the irrigated area of the mth municipality in 1999 and 2009 (in ha), respectively and IA RA99 is the irrigated area of its agrarian region (in ha).The extrapolated irrigated areas for 1999 were further disaggregated into different crops and irrigation systems by using the respective fractions of 2009 compared to the total irrigated area of the same year.Crop N application rates from sludge, manure and inorganic fertilisers were derived from Serra et al. (2019) and distributed to crop areas under different irrigation systems.
The MITERRA-EUROPE approach (Velthof et al., 2009) was used to estimate runoff and leaching N losses (Table 1).Data regarding the GNB and gaseous losses were spatially distributed to a spatial resolution of 100 × 100 m using the agricultural areas derived from the Corine Land Cover (CLC) closest year, 2000 and 2012.An adjustment factor was used to match the acreage of the utilized agricultural areas at the municipality scale from Statistics Portugal (1999Portugal ( , 2009) ) and CLC agricultural areas for the nearest year.N losses in surface runoff were computed as: where N app is the N-input from sludge, manure spreading, inorganic fertilisers and animal N excreted onto pastures, corrected to NH 3 losses during housing, storage and grazing (kg N ha −1 yr −1 ) and f runoff is the runoff fraction (%N app ).The latter is obtained as a function of slope (S) and land use (LU) modified according to the excess of precipitation surplus (PS; precipitation minus evapotranspiration), soil type (T) and depth to rock (DR), reclassified according to Velthof et al. (2009): NO x emissions from soils EMEP (2016) Tier 1 Application of N in manure, sludge, fertiliser and excreted onto pastures (Serra et al., 2019).
N leach was quantified as the product of the SSNB and leaching fractions (f leach ; %SSNB): where f leach is determined as a function of soil texture (ST) and LU modified according to PS, soil organic carbon (SOC), temperature (TE) and rooting depth (RD) (Eq.( 7); Velthof et al., 2009).The remaining fraction of the SSNB (f denit ; %SSNB) concerns denitrification to N 2 .

Estimating groundwater N contamination from diffuse agricultural losses
In mainland Portugal, there are 93 groundwater bodies aggregated into four main hydrogeological units (HU): Hercynean massif, Meridional, Western and Tagus-Sado (Fig. 1).Sixty are aquifers, herein mentioned as the main aquifer systems, and 33 are undifferentiated hydrogeological formations.The main aquifer systems are located in the coastal regions of Central (Western unit and Tagus-Sado) and Southern Portugal (Meridional unit) (Fig. 1).The Hercynean massif, mainly composed by igneous and metamorphic rocks, is characterized by low permeability and productivity except in some carbonate aquifers in the Alentejo province (Ribeiro and Cunha, 2010).The Meridional unit contains carbonate aquifers that are slightly over-exploited for agricultural purposes.The Western unit has sedimentary and karst aquifers as well as some multi-layered aquifer systems.The Tagus-Sado basin is the most important aquifer in the Iberian Peninsula and the main source of water for societal and agricultural purposes (Almeida et al., 2000).It is a multi-layered aquifer system, covering 8000 km 2 in a sedimentary basin where intensive agriculture occurs with implications in groundwater pollution (Cameira et al., 2019(Cameira et al., , 2021)).

Water percolation
Aquifer recharge results from percolated water from precipitation and irrigation.Aquifer recharge rates presented in the River Basin Management Plans for the eight hydrographic regions in mainland Portugal were derived only from precipitation (RR; in % precipitation; Lobo-Ferreira et al., 2015).Since irrigation considerably impacts aquifer recharge in Mediterranean regions, especially during the summer, here the recharge was estimated from both precipitation and irrigation data.Thus, the water that percolates below the root zone (Perc; m 3 ha −1 yr −1 ) is used as a proxy of aquifer recharge, similarly to other studies (Arauzo, 2017;Cameira et al., 2021).It was calculated as the result of the water balance equation (Eq.( 8)) where P is the precipitation (mm yr −1 ), I the irrigation water (m 3 ha −1 yr −1 ), ET c the crop evapotranspiration (m 3 ha −1 yr −1 ) and R is the water lost through runoff (m 3 ha −1 yr −1 ).The irrigation water supply (I i,j ; m 3 ha −1 yr −1 ) to the ith crop irrigated by the jth system and occupying the area IA i,j (in ha) at the municipality scale was estimated based on (i) historical regional crop irrigation requirements (IR i,j ; m 3 ha −1 yr −1 ) collated from the official national irrigation authority (DGADR, 2019) and (ii) the irrigated areas per crop and irrigation system previously used in the direct N 2 O calculations.Data regarding crop acreage was spatially distributed according to different land uses (LU) from the Corine Land Cover database (Permanently Irrigated, Rainfed, Olive groves, Orchards, Pastures, Rice, Vineyards, Heterogeneous).The total irrigation amounts for the ith specific crop (I i ; in m 3 yr −1 ) was computed as: Crop evapotranspiration was calculated according to the FAO crop coefficient method (Allen et al., 1998): where ETo (m 3 ha −1 yr −1 ) is the reference evapotranspiration and K c,LU is a crop-specific coefficient variable according to the growth stage (dimensionless), aggregated for each LU class.Finally, runoff losses were calculated as a function of precipitation, irrigation, crop evapotranspiration and runoff fraction (f runoff ; v.d 2.2): 2.3.2.Risk of aquifer contamination with nitrate from agriculture An index-based method was used to quantify the potential risk of groundwater contamination by agricultural diffuse N losses.This risk index (RI; dimensionless) is obtained by overlaying the N contamination hazard with the aquifer vulnerability (Cameira et al., 2021;Kazakis and Voudouris, 2015).The N hazard, represented by the nitrate concentration in the potential aquifer recharge from below the root zone (N c ; mg NO 3 − L −1 ), improves on previous approaches that used the N surplus (Cameira et al., 2019) and land use ratings (Kazakis and Voudouris, 2015).Aquifer vulnerability is represented by residence time (RT; yr), soil texture and slope, incorporating the effect of intrinsic vulnerability (groundwater depth, vadose zone lithology) and specific vulnerability (aquifer recharge).Nitrate concentration in the leached water (N c ) was calculated as the ratio of the total N-leaching (N leach ; kg N yr −1 ) from Eq. ( 6), and the percolation below the root zone (Perc; m 3 yr −1 ) from Eq. ( 8).Nitrate residence time in the vadose zone was calculated according to Ascott et al. (2017): where D is the groundwater depth (m), V NO3 is the NO 3 − vertical average velocity from the bottom of the rootzone to the groundwater (m yr −1 ) and ϕ is the effective porosity of the vadose zone (dimensionless).
Groundwater depth and effective porosity were collated from Fan et al. (2013) and Gleeson et al. (2011), respectively.Slope and soil texture were obtained from the same sources as described in Table 1.All parameters were resampled to a spatial resolution of 100 × 100 m and reclassified into five scores ranging from <0.2 (low risk) to 1 (very high) (Table 2).Vulnerability was computed as a function of residence, slope and soil texture with three different relative weights for each vulnerability parameter: where w rt , w slope and w soil texture are the weights for residence time (RT), slope and soil texture, respectively, ranging from 0 to 1.The relative weights for each parameter were computed emulating the approach to calculate "effective weights" in single parameter sensitivity analysis.This approach is widely used to assess the influence of the vulnerability parameters of the DRASTIC method (Krogulec and Trzeciak, 2017;Oke, 2020;Rajput et al., 2019;Tomer et al., 2019).Firstly, vulnerability was calculated by setting w rt to 0.6, and w slope and w soil texture to 0.2, as the theoretical weights.A higher weight was given to RT since the N leached already accounts for soil texture and slope from the MITERRA fraction parametrization (i.e., f leach ).Moreover, the time lag until nitrate seeps into the phreatic zone depends on residence time (RT), representing a hydrologic legacy effect (Chen et al., 2018).Moreover, the time lag until nitrate seeps into the phreatic zone depends on residence time (RT), representing a hydrologic legacy effect (Chen et al., 2018).Spatially explicit effective weights for the ith vulnerability parameter (W i ) were then estimated: where P w and P r are the theoretical weight and reclassified raster map for the ith parameter and V is the vulnerability.All units are dimensionless.The vulnerability index was then recalculated using the effective weights.The first iteration of vulnerability and the effective weights are displayed in the supplementary material (Figs.S1 and S2).
The Risk Index (RI) for groundwater contamination with agricultural nitrates can be thus computed: RI was assessed following the risk score presented in Table 2.All (GIS-based) calculations were performed using a spatial resolution of 100 × 100 m using R version 3.6.2(R Core Team, 2019).

Overview of the soil surface nitrogen balance in Portugal
There was a small reduction of gaseous and runoff losses between 1999 and 2009 (Fig. 2, left).Volatilization of NH 3 was the predominant pathway of N loss, averaging 30 kg N ha −1 in both years, with a slight decrease in maximum values (292 to 282 kg N ha −1 yr −1 ).Emissions of other gases (N 2 O and NO x ) were relatively small in most municipalities, with average values below 2.5 kg N ha −1 yr −1 .The direct N 2 O emissions were approximately 15% lower when using the Mediterranean EFs (0.75%) comparatively to the IPCC EFs for temperature regions (1%).Despite the similar runoff fractions in both years (≈4.8%),runoff losses decreased on average almost 45% (3.8 to 2.1 kg N ha −1 yr −1 ).
Fig. 2 (right) shows the dispersion of the SSNB calculated for all municipalities.The overall value for the country decreased 26% from 1999 to 2009 (123 to ~91 Gg N yr −1 , respectively) (1 Gg = 1000 t).However, the variability among municipalities remained large and even increased (more outliers), nearly reaching 600 kg N ha −1 yr −1 in some municipalities in the coastal areas of the Central and NW regions.

Estimating nitrogen leaching from the soil surface nitrogen balance
Table 3 presents the leaching and denitrification fractions, as percentage of the soil surface N balance (SSNB) and the corresponding amounts of N leached from below the root zone and through denitrification to the atmosphere at the mainland level.The largest portion of the SSNB was allocated to denitrification (≈81%).Total N leached decreased 32% over time, from 23.5 to 15.9 Gg N yr −1 .This occurred in spite of setting N-leaching to 0 in 23 municipalities in 1999 and 9 in 2009 due to negative SSNBs, corresponding to 8.3 and 3% of the 278 municipalities, respectively.These were mainly located in NW Portugal (Fig. 3), some within the upper section of the Western unit, and totaled 244 t N leached in 2009.Fig. 3 shows the spatial variation of N-leaching from below the root zone calculated at a grid of 100 × 100 m in 1999 and 2009.On average, from 1999 to 2009, N-leaching decreased from 5.5 to 3.7 kg N ha −1 yr −1 , particularly in the South, over most of the west coast and to a lesser extent in the inland regions.By contrast, some coastal areas in central Portugal showed higher annual N-leaching losses (50.7 to 68.0 kg N ha −1 ) in 2009 compared to 1999.Despite the overall decrease in N leaching in mainland Portugal, some parts of the Tagus-Sado basin and Western unit showed an increase in the N leached (Fig. 3).

Percolation water
Percolation water below the root zone is explored as an estimation of the recharge.In agricultural land, the potential aquifer recharge (Perc) includes both the precipitation and irrigation water.At the mainland level, Perc approximately halved between 1999 and 2009, declining from 8.8 to 4.4 km 3 yr −1 .On average, Perc also showed an average reduction of 50%, from 2067 ha −1 yr −1 to 1043 m 3 ha −1 yr −1 , respectively.Perc declined over most of mainland Portugal but particularly in southern Portugal but also to a lesser extent in some regions of NW Portugal (Fig. 4).These reductions occurred despite an increase of 21% in irrigation water usage (1.8 to 2.2 km 3 yr −1 ), mainly in west coast of Portugal and the Alentejo region.Fig. 4 also shows the variation between 1999 and 2009 of Perc aggregated at the aquifer level.Perc in the main aquifer systems totaled 1.8 and 0.9 km 3 for 1999 and 2009, respectively.This decline  resulted in lower estimates for 49 of the 60 aquifers in mainland Portugal with an average reduction of 70%, while the remaining aquifers showed comparatively lower increases (20%).The reduction in Perc was greatest to the west and south of the country, with the increases occurring mainly to the north of this area, near the coast.

Risk assessment of groundwater N contamination or pollution from diffuse sources
Fig. 5 displays the spatial variability of the hazard, vulnerability and the risk index in 1999 and 2009.On average, the hazard increased from a moderate to a high level (mean of 0.41 in 1999 and 0.69 in 2009), with the area associated with a very high hazard increasing from 0.17 to 0.63 million ha over the ten years.There was a considerable change in the spatial distribution of cells with a very high hazard.In 1999, the high/ very high hazard areas were scattered throughout the Hercynean massif aquifers in the Alentejo with some minor hotspots in the other HUs.In the Tagus-Sado basin, the hazard was moderately high in its upper section with some cells achieving a high hazard while its lower section achieved a lower hazard.However, in 2009 most of this HU attained a very high hazard.There was an increasing tendency in southern Portugal, while the upper part of the Western unit remained similar in both years.
Vulnerability was composed of one temporal (residence time) and two physical (slope and soil texture) parameters.Both parameters were higher in the sandy soil flatlands of southern Portugal and the Tagus-Sado Basin (Fig. S3).In the upper section of the Western unit the slope parameter was very high, but the soil texture parameter was low, while the opposite occurred in its lower section.Residence time increased from 1999 to 2009 (Fig. S4), except in the upper section of the Western unit whose partial vulnerability remained very high in both years.On average, the effective weights (i.e., contribution to the vulnerability) were greater for residence time (50 and 41%), slope (30 and 34%) and soil texture (21 and 25%) for 1999 and 2009, respectively.
The total vulnerability was, on average, classified as high and remained very similar in both years although with a slight declining tendency (0.69 in 1999 and 0.64 in 2009).There was a reduction in the areas classified as very high (0.24 to 0.19 million ha) but with a small increase in areas with a high vulnerability (0.36 to 0.39 million ha).Vulnerability was higher in the upper and main sections of the Western unit and Tagus-Sado, respectively, and in clusters throughout southern Portugal.
The groundwater contamination risk index (RI) increased, on average, from a Low (0.29) to a Moderate risk (0.44).There was a considerable expansion of areas with very high risk, which grew from 8800 ha in 1999 to 82,679 ha in 2009.This increase occurred mainly in southern Portugal, particularly in the Tagus-Sado basin, the central and right side of the Meridional unit and in all aquifers of the Hercynean massif.The aquifers with the highest area fraction with very high risk were located in the Hercynean massif: the Moura-Ficalho aquifer (20% or 3401 ha) in 1999 and the Gabros de Beja aquifer (40% or 12,589 ha) in 2009.The three aquifers of the Tagus-Sado basin attained the largest areas classified as very high risk in 2009 (>28,000 ha each).Conversely, the spatial variation in the Western unit remained similar and relatively low in both years, especially as most of its upper section was classified with a very low risk.There was also an increase from a very low to a very high risk in a minor hotspot in its southern section.

Gaseous and runoff losses
To gain an insight into the spatial distribution of the main N losses related to agricultural activities requires an understanding of the dynamics of land cover and land use in mainland Portugal.Between 1999 and 2009, the total arable land decreased by 30% (1.7 to 1.2 million ha) whereas the acreage of extensive pastures increased by a similar fraction (30%; 1.3 to a 1.7 million ha) (Statistics Portugal, 2009).This is partly explained by contrasting land-use tendencies between inland and coastal regions.The former are characterized by land abandonment and the conversion of arable land to extensively managed pastures and agroforestry systems (Malek and Verburg, 2018).In contrast, the coastal regions, where the most productive hydrogeological units are located (Fig. 4), experienced an improvement in infrastructures (e.g., roads, hydraulic infrastructures), which led to the expansion of peri-urban agriculture (maize and vegetables) and livestock production (poultry and pigs).The concentration of agricultural inputs and livestock into progressively smaller areas have had a direct impact on nitrogen pollution and water quality.
Although gaseous/runoff losses decreased at the mainland level, the number and magnitude of the outliers increased (Fig. 2).The increase of local N losses is particularly evident in the coastal regions of Central Portugal (upper section of the Western unit).The main drivers for this were the increase in the rates of manure spreading (393 to 600 kg N ha −1 yr −1 ), which increased runoff losses, and an increase in urea-based fertilisers (14.5 to 24.1 Gg N yr −1 ) that led to local higher NH 3 emissions, despite the national reduction in inorganic fertiliser N (149 to 96 Gg N yr −1 ; APA, 2017).These areas were identified as agricultural N pollution hotspots by Serra et al. (2019), who found GNBs up to 1000 kg N ha −1 yr −1 despite the overall decreasing GNB tendency at the mainland level, from 49 to 39 kg N ha −1 yr −1 between 1999 and 2009.The similar pattern was found in SSNB which suggests a higher potential for local losses of N leaching.

Nitrogen leaching
The reduction of N leach in Portugal between 1999 and 2009 (32%) derived from the overall lower leaching potential as estimated by the SSNB and from lower f leach (Table 2), mostly due to the decline in precipitation.Both f leach (18%SSNB) and f runoff (5%N-inputs) are in line with estimates from Keuskamp et al. (2012) (28 and 5%, respectively) and Velthof et al. (2009) (3 and 12%, respectively) for Europe.Similarly, our mean value for N leach in Portugal for 1999 is quite comparable to data from Velthof et al. (2009) for the year 2000 (6 and 8 kg N ha −1 yr −1 , respectively).The decrease in UAA may have contributed to this reduction of N leach in mainland Portugal, but other factors are also likely to be important.These include the implementation of the ND and the delineation of Nitrate Vulnerable Zones (NVZ) in 2004 through Action Programs, with mandatory measures such as restrictions on animal stocking rates and maximum fertilization rates per crop and soil physical properties.According to Velthof et al. (2014), the implementation of the ND in Europe reduced N leaching losses by 16% for the period 2000-2008, but up to 30-60% in countries with intensive agriculture such as in northwestern Europe where the NVZs cover larger areas (Dalgaard et al., 2014;Van Grinsven et al., 2012).
In mainland southern Portugal there was a considerable improvement in reducing N leach , which may be explained by land abandonment/conversion to extensively managed permanent pastures, a greater efficiency in the use of inorganic fertiliser N, and/or lower inputs of manure N although it is unclear which is the predominant factor (Serra et al., 2019).Our estimates for one NVZ in this region (240 t N in 2009) are comparable to those also made by Malta et al. (2017) and Hugman et al. (2017) (300 and256 t N for 2007, respectively).By contrast, N leach increased in the remaining aquifer systems, particularly in the upper section of the Western unit and in the Tagus-Sado basin to a lesser extent (Fig. 3).However, N leach was set to 0 in some regions within the limits of two NVZs located in the Western unit.Cameira et al. (2019Cameira et al. ( , 2021) ) and Cordovil et al. (2018) present similar conclusions for the spatial distribution of the Tagus NVZ insofar as there was a significant reduction of the GNB in its upper section and a slow stabilization in its southern part, where livestock production is the main activity.Additionally, the expansion of intensively managed irrigated agriculture (41,962 to 48,113 ha) may have contributed to this diverging spatial tendency.
The implementation of the NVZs, which cover only 4.4% of mainland Portugal, resulted in different levels of effectiveness in reducing N leach .Our data suggests that regions with low manure N (Algarve in the Meridional unit; Fig. 1) or where grazing practices are prevalent (Alentejo in the Hercynean massif; Fig. 1) were able to achieve greater reductions in N leach .Conversely, the Western unit and the Tagus-Sado basin, where there is a greater competition for land, resulted in the intensification of local agro-ecosystems and subsequently of aggravated nitrate pollution.This occurred despite the implementation of three NVZs, two in the Western unit and one in the Tagus-Sado basin.

Percolation water
The year 2009 was comparatively much drier than both the mean historical data and 1999 (Fig. 1), registering a reduction of 50% in Perc over mainland Portugal and also in the area covered by the main aquifer systems.Accordingly, our estimates show that irrigation water requirements increased over time (1.8 to 2.2 km 3 yr −1 ) and are in line with data from FAO (2016), around 2 km 3 for the year 2007.Irrigation was of particular importance in the Western unit and Tagus-Sado, partly due to the existence of rice paddies and the intensification of agricultural production (e.g., tomatoes for industrial processing and forage maize), but also in the Alentejo region in 2009, also as a result of the access to irrigation water from the new Alqueva reservoir.However, irrigation could not counteract the comparatively lower precipitation in 2009, resulting in a sharp decline in Perc, particularly in southern Portugal (Fig. 4).The marked spatial variation of Perc in the upper section of the Western unit derived from a small increase of 10 and 25% in precipitation and irrigation, respectively.

Identification of areas with high risk of agricultural nitrate contamination
Our results show a large increase in agricultural areas identified as having a very high hazard and groundwater contamination risk despite the overall reduction in nitrate leaching from below the root zone (32%).This increase derived from a larger decline in Perc (50%) which resulted in the overall increasing tendency for nitrate concentration in percolation water from below the root zone (i.e., hazard).Locally, the increase in hazard holds for two different situations: one where both N leach and Perc increased, albeit where the latter did so at a faster rate (e.g., Southern Portugal), and other where N leach increased and Perc decreased (e.g., Tagus-Sado basin).Conversely, the upper section of the Western unit showed an increasing tendency for N leach and Perc, thus resulting in a slightly higher hazard.It should be noted that the hazard in the upper section of the Western unit may have been underestimated, since N leach was set to 0 in some municipalities, due to negative GNBs (Serra et al., 2019).
This approach to estimate hazard fails to account, however, for nitrate accumulation in the vadose zone from past agricultural activities (Ascott et al., 2017).While the overall reduction in Perc increased the agricultural hazard, it had an opposite effect in RT, which was the vulnerability parameter with the highest effective weights in both years.The Western unit attained a very high vulnerability in both years as most of the area is characterized by flat, sandy soils with low residence times.Our estimates of areas with very high vulnerability are in agreement with studies in the Western unit (Ribeiro and Cunha, 2010) and the alluvial plains in the Tagus-Sado basin (Mendes and Ribeiro, 2010) using a Susceptibility Index.
The diverging tendencies calculated for the hazard and vulnerability resulted, nonetheless, in an increase of the overlap between cells with high hazard and vulnerability, and thus in the risk of nitrate contamination (Fig. 5).The reduction in Perc, particularly in southern Portugal, may have over-and underestimated the hazard and vulnerability, respectively.The areas classified with very high risk in 2009 closely resembled the NVZs in Portugal, such as the Tagus NVZ in the Tagus-Sado basin, the Gabros de Beja NVZ in the Hercynean massif aquifers or the Tavira NVZ in the eastern part of the Meridional unit (Fig. S5).Furthermore, the development of this risk index allowed the identification of a new hotspot in the central region of the Meridional unit.This region is characterized by intensively managed citrus and vegetable productions (6156-8730 m 3 ha −1 yr −1 for drip irrigation, respectively; DGADR, 2019) and with inputs of inorganic fertiliser N (ca.130 kg N ha −1 yr −1 ).Despite the increase in the contamination risk, our data suggests an improvement in N leach (Fig. 3) and in the application of inorganic fertilisers (Serra et al., 2019) in this region.It should be noted, however, that by not accounting nitrate storage in the vadose zone, the contamination index fails to identify the risk derived from the cumulative effect of relatively low hazard but very high vulnerability areas (e.g., upper section of the Western unit).Further, the index may overestimate the risk in regions with high hazard but a very low vulnerability since nitrate transport from the groundwater to aquifers is reduced due to the intrinsic conditions.Therefore, these contrasting cases may require different targeted measures to reduce the contamination risk (see below).

Usefulness in the definition of spatially targeted measures to tackle groundwater contamination
The risk index here developed provides a relatively simple approach that enables the identification of groundwater vulnerable areas, such as by other indexes (e.g., DRASTIC; USEPA, 1985), while its scope is extended by including a N hazard module.For this reason, it presents an alternative approach to delimitate NVZs (e.g., Arauzo, 2017;Sajedi-Hosseini et al., 2018), and an alternative to process-based models, which require considerably more input data.Its application improves our understanding of the interactions between hazard and vulnerability, integrating these into a framework to assess the risk of groundwater contamination.Furthermore, it contributes to groundwater governance by simplifying information flows about prevalent issues and thus to enact countermeasures accordingly.
The index is particularly useful when applied to contrasting situations of vulnerability and hazard, which require distinctly different measures.Areas with high vulnerability but low hazard (e.g., upper section of the Western unit and Tagus-Sado basin) may require drastic interventions to preserve the quality of groundwater resources, since N leached easily reaches the phreatic zone.Spatially differentiated strategies regarding changes in land use and land cover can be employed to reduce the hazard (see Hashemi et al., 2016 for a comprehensive review).For instance, conversion of cropland to pasture or forest conversions can reduce N leaching (Li et al., 2020;Liu et al., 2013;Rode et al., 2009) and enhance soil water storage (Kroes et al., 2019), thus positively affecting both the hazard and vulnerability.However, these measures may have limited effectiveness when employed in areas with high vulnerability.This is illustrated by the Netherlands, a country with relatively high groundwater vulnerability (Nistor, 2019), where despite increasingly stringent nutrient policies, the groundwater nitrate concentrations still exceeded 50 mg L −1 in the sandy regions during the period 1990-2014(Van Grinsven et al., 2016).In these regions, reducing the groundwater contamination risk and nitrate concentration may require such severe restrictions on agricultural activities that they become untenable.If this is the case, the reallocation of agricultural activities to areas with lower vulnerability would be an option.However, this would have severe socioeconomic implications and it is unlikely to be implemented unless top priority is given to groundwater governance.
In contrast to areas with a high vulnerability and low hazard, areas with low vulnerability and high hazard (e.g., southern Portugal) can greatly benefit from stricter agri-environmental policies and associated agronomic measures, through the expansion of the designated NVZs.However, our data shows that the positive effect of the implementation of some existing NVZs did not achieve the expected effect, due to the reduction in precipitation.Improved agronomic practices in Mediterranean agro-ecosystems can reduce seasonal nitrate leaching losses through cover cropping during autumn/winter heavy precipitation events (Poch-Massegú et al., 2014) and optimized N fertilization and irrigation practices in the summer (Sanz-Cobena et al., 2017).The use of nitrate-contaminated groundwater for irrigation presents an additional source of N, often underestimated by farmers.In some parts of the Tagus-Sado basin, irrigation N accounted for 30-35% of total N-inputs (Cameira et al., 2019), though there are reported values up to 159 kg N ha −1 yr −1 elsewhere (Quemada et al., 2020).Shifting the irrigation source to surface water had a greater impact in nitrate concentrations than improved agronomic practices in an aquifer in the Meridional unit (Hugman et al., 2017), although its effectiveness as a general measure would vary, depending on nitrate concentration in surface-and groundwater.

Future considerations
In the future, the model here employed can be expanded over an annual basis for a larger period to better understand temporal trends, with the inclusion of the Portuguese islands (Azores and Madeira), which have a different set of hydrogeologic conditions and N dynamics (Prada et al., 2016;Cruz et al., 2013).Other anthropogenic point sources (e.g., sewage effluents, industry) that affect groundwater quality can be calculated (e.g., Morée et al., 2013) or derived from public datasets (Vigiak et al., 2020).Similarly, the flexibility and scope of the vulnerability modelling framework can be enhanced by including other nutrients and elements than nitrogen.For instance, Kazakis et al. (2020) and Busico et al. (2020) collected water and soil samples for multiple parameters used as vulnerability predictors.Although Busico et al. (2020) achieved moderately strong correlations with nitrate concentrations (R 2 of 0.62 and 0.75), machine learning algorithms such as random forests also have the potential to accurately predict nitrate concentration in groundwater (Canion et al., 2019;Knoll et al., 2019;Rahmati et al., 2019).A scenario module can also be developed to gain insight on future trends of agricultural N flows and associated impacts in terms of risk areas (Nistor, 2019).This can provide a valuable contribution to (ground)water governance.

Conclusions
The risk index of nitrate contamination proposed here aimed to assess the groundwater contamination risk with N from agricultural activity, at a national scale, considering the spatial variation of the agricultural N hazard and the site-specific aquifer vulnerability.As an input to the contamination risk, percolation has opposite effects on the hazard and the vulnerability: its reduction can increase the agricultural hazard while reducing the vulnerability.

Fig. 2 .
Fig. 2. Agricultural gaseous (NH 3 , N 2 O, NO x ) and runoff N losses (left) and the soil surface N balance, SSNB, the proxy to estimate N-leaching (right).The left plot uses a logarithmic scale.

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3. Total N leaching from below the root zone for 1999 and 2009, and the spatial variation between 1999 and 2009 in the main aquifer systems.

Table 1
Overview of the different data requirements and methods employed in the estimates of the various environmental N losses.

Table 2
Risk scores used in the reclassification of N c and RT

Table 3
Denitrification and leaching fractions of the Soil Surface N balance (SSNB), and total N denitrified and leached in mainland Portugal in 1999 and 2009.The numbers in brackets are given in kg N ha −1 yr −1 .