Assessing Natural Background Levels in the Groundwater Bodies of the Apulia Region (Southern Italy)

: Deﬁning natural background levels (NBL) of geochemical parameters in groundwater is a key element for establishing threshold values and assessing the environmental state of groundwater bodies (GWBs). In the Apulia region (Italy), carbonate sequences and clastic sediments host the 29 regional GWBs. In this study, we applied the Italian guidelines for the assessment of the NBLs, implementing the EU Water Framework Directive, in a south-European region characterized by the typical Mediterranean climatic and hydrologic features. Inorganic compounds were analyzed at GWB scale using groundwater quality data measured half-yearly from 1995 to 2018 in the regional groundwater monitoring network (341 wells and 20 springs). Nitrates, chloride, sulfate, boron, iron, manganese and sporadically ﬂuorides, boron, selenium, arsenic, exceed the national standards, likely due to salt contamination along the coast, agricultural practices or natural reasons. Monitoring sites impacted by evident anthropic activities were excluded from the dataset prior to NBL calculation using a web-based software tool implemented to automate the procedure. The NBLs resulted larger than the law limits for iron, manganese, chlorides, and sulfates. This methodology is suitable to be applied in Mediterranean coastal areas with high anthropic impact and overexploitation of groundwater for agricultural needs. The NBL deﬁnition can be considered one of the pillars for sustainable and long-term groundwater management by tracing a clear boundary between natural and anthropic impacts.


Introduction
The EU Water Directives [1][2][3], implemented in Italy by Acts [4][5][6], require Member States that good status of groundwater bodies (GWBs) is achieved or maintained to protect all dependent ecosystems. Chemical status is determined based on quality standards set in Europe for nitrates and pesticides, and threshold values (TVs) that have to be set by Member States for all contaminants or groups of contaminants characterizing the GWBs as being at risk of failing the environmental objectives. The NBL can be defined as the concentration of a substance, or the value of an indicator, in a given GWB or group of GWBs, with no or extremely limited anthropogenic modifications. In this perspective, the TV of a considered natural substance may be redefined by the local authorities based on the NBL. The EU Directives do not provide any methodology for the TVs and NBLs establishment in the groundwater and the Member States are left free to apply their own approach considering the conceptual models of the groundwater bodies.
Chemical composition of groundwater depends on a wide range of natural factors such as water-rock interaction, residence time, rainfall input, chemical and biological processes that occur in the unsaturated and saturated zone, interactions with other water bodies [7]. Distinguishing between natural and anthropogenic sources is not a simple task since it requires high-level knowledge of all processes involved and the availability of a proper amount of monitoring data. The simplest methodologies for the NBL estimation within a GWB are based on the temporal and spatial approaches or on a combination of them [8,9]. The temporal approach uses historical groundwater quality data, which are assumed to be representative of the existing conditions before supposed human impacts. However, historical chemical data are usually scarce and have been determined using different methods and protocols compared to the modern sampling and analytical standards. The spatial approach uses monitoring data collected from areas in which the anthropogenic activity is assumed low and with no impact on the water quality. In this case, some care should be paid, when extending the NBLs values calculated at unspoiled areas to those impacted since they could have different physical, chemical, and biological characteristics. The difficulty of finding historical datasets and identifying pristine portions of the GWBs in populated areas has often led to developing other approaches.
The guidance issued within the EU 6th FP BRIDGE project [10,11] suggests methods based on either a statistical or a pre-selection approach. The partition of populations method [12][13][14][15][16] is a statistical approach based upon the idea that the concentration distribution of a chemical species in a groundwater system can be considered as the combination of two or more components of natural and anthropogenic origin.
On the other hand, a pre-selection based approach eliminates samples clearly affected by human activities on the base of specific markers such as nitrate, ammonium, chloride, synthetic organic compounds, etc. The samples exceeding a fixed value of the chosen markers are eliminated and the NBL is chosen as the 90th, 95th, or 97.7th percentile of the remaining water samples dataset. The main disadvantage of this approach is that the samples are eliminated according to qualitative criteria [17], e. g. the amount of available data and the shape of the related statistical distribution, and that the removal of samples might result in a too scarce data set.
Urresti-Estala et al. [18] assessed the NBLs using the objective criteria of two statistical techniques: the iterative 2σ technique, and the distribution function. Differently from the commonly employed statistical techniques, these Authors proposed two different approaches that do not require the elimination of samples, nor specific prior distributions, or the presence of large data sets. Both techniques are capable to discriminate between high values due to natural background populations from those induced by human influence.
Different attempts to apply a multi-method approach for the NBLs assessment have been made generally combining the pre-selection method and statistical techniques [19][20][21][22][23]. These studies pointed out that their integration reinforces the validity of the assessment, particularly when the hydrogeological setup and geochemical characteristic of groundwater are properly considered, confirming the greater importance of the conceptual model of the system than the choice of a statistical method. Molinari et al. [24] proposed other methodologies integrating partition population and numerical modeling of flow and transport processes to estimate NBLs of manganese and sulfate in potentially contaminated coastal aquifers. This study highlighted the need to validate the NBLs assessment procedures with accurate hydro-geochemical modeling results to identify external forcing components influencing the NBLs, such as the effects of tides and seawater intrusion, and avoid the possibility of interpreting naturally induced concentrations as anthropogenic effects.
Recently, IRSA-CNR collaborated with the Italian Institute for Environmental Protection and Research (ISPRA) in defining scientifically based guidelines for the NBLs assessment and clarifying some methodological aspects [25]. These guidelines rely on the pre-selection method but leave to the local water managers the task of identifying the specific pressure indicators while the only mentioned parameters are nitrate (for oxidizing conditions) and ammonium (for the reducing environments). Further, the guidelines do not mention how to deal with saline groundwater; hence, a procedure for the definition of a boundary between contaminated and uncontaminated groundwater in coastal aquifers has still to be defined. Assessing the origin of groundwater salinization in coastal aquifers is an issue mainly because it can depend on many concurrent natural or anthropic causes. Besides the possible contribution of seawater intrusion driven by overexploitation and sea-level rise, the scientific literature reports a detailed list of potential sources [26,27].
This work proposes a tentative application of the national guidelines for the assessment of the NBL values of some inorganic compounds in the GWBs of the Apulia region (southeast of Italy). In the case study, nitrate and ammonium concentrations have been considered as marker of anthropic contamination of the water samples. This work has been carried out within the framework of the research project on the assessment of background levels in the groundwater bodies of the Apulia region (VIOLA) jointly developed with the Department of Water Resources Management of the Apulia Region.

Materials and Methods
The Italian guidelines for the NBLs assessment describe an integrated approach based on the pre-selection and statistical methods [25]. The procedure consists of ten steps, some of which (steps 1-4) apply on the monitoring sites (MSs) and others (steps 5-10) on the single parameters: 1.
Conceptual model definition.

2.
Data exploration (validation of analytical data and handling of non-detect data).

3.
Hydrochemical facies identification, based e.g., on dissolved oxygen concentration and/or redox potential. The facies analysis, when performed, leads to the separation of the dataset into separate subsets (e.g., reducing/oxidizing datasets) with different geochemical characteristics to be processed separately.

4.
Pre-selection of data using the chosen anthropic markers (e.g., NO − 3 concentrations in oxidizing facies and NH + 4 concentrations in reducing facies); the resulting data are compiled in the so-called "pre-selected dataset". 5.
Analysis of the time series per monitoring site (MS): identification and handling of the outliers by graphical (box and whisker plot) and statistical method (Huber's test), trend analysis, selection of the median (or single value when only one measure is available) to be used as "representative value" in the following step. 6.
Analysis of the data at the GWB scale (one data per monitoring station): identification and handling of the outliers of the preselected datasets, assessment of the existence of multiple populations and consequent identification of sub-sets. 7.
Analysis of the statistical distribution (using the Shapiro-Wilk test) of the target parameter(s). 8.
Assessment of pre-selected dataset size (significant spatial size: monitoring stations ≥ 15, significant temporal size: at least 8 half yearly observations); 9.
Identification of specific path types based on pre-selected dataset size; 10. Assessment of the NBLs by statistical methods according to the specific path.
The pre-selected dataset size is a key element because it determines the procedure for NBL assessment and the robustness of the results. The guidelines identify four possible scenarios of dataset size (see point 9 of the previous list) considering both spatial and temporal dimension ( Table 1).
The obtained NBLs are characterized by different confidence level (high, medium, low, very low) with regard to: the number of total observations or the number of total monitoring sites, the extension of GWB, and the aquifer type (confined or unconfined) ( Table 2). When the confidence level is low or very low, the NBL is considered provisional. The confidence level is used to select a NBL at GWB scale from those calculated for different facies.  The logical phases and the requirements of the guidelines have been implemented at the IRSA-CNR into an automated tool capable of standardizing NBLs assessment operations. This tool, eNaBLe [28], thanks to its modularity, has been successfully integrated into the VIOLA project web site (http://www.irsa.cnr.it/Viola/, accessed on 30 March 2021) which hosts different statistical, graphical elaboration, and visualization procedures. eN-aBLe and the other VIOLA web site procedure are written in PHP [29] and use MySQL [30] as relational databases manager.
The database is structured in 5 main tables (Monitoring Points, Wells, Springs, Samples, Chemical Analyses), 2 ancillary tables (GWBs, Parameters) and one table for user authentication and privileges. The user interface of the eNaBLe tool is organized in three logical phases corresponding to three different associated web pages: (i) Selection of the GWB data to be processed, analytical parameters to take into account and selection of the calculation options to be used; (ii) NBLs assessment and relative confidence levels; (iii) output of results in graphical, tabular and report form.
By performing a query on the selected dataset, eNaBLe produces a preliminary table in which all the analytical parameters and their basic statistical description are listed, proposing for the NBLs evaluation those analytical parameters in which at least one exceedance of the relative TV has been found. Next, it is possible to select the options for the NBL assessment, such as the appropriate time interval for the data to be processed, the parameters to be used for the redox facies separation and the limits to apply for the preselection process to nitrates or ammonium in relation to the oxidizing or reducing facies, among others.
A screenshot of the GWB selection and parameters configuration page is shown in Figure 1.
By performing a query on the selected dataset, eNaBLe produces a preliminary table in which all the analytical parameters and their basic statistical description are listed proposing for the NBLs evaluation those analytical parameters in which at least one exceedance of the relative TV has been found. Next, it is possible to select the options for the NBL assessment, such as the appropriate time interval for the data to be processed the parameters to be used for the redox facies separation and the limits to apply for the preselection process to nitrates or ammonium in relation to the oxidizing or reducing facies, among others.
A screenshot of the GWB selection and parameters configuration page is shown in Figure 1. The results of the calculation as well as the settings used are listed in the second page of the procedure (Figure 2). In addition, the result page provides the links to the intermediate products of the calculation such as the validated dataset (built after time filtering, non-detected values handling and electrical balance verification), pre-selected and facies separated datasets, final datasets used for the calculation and the path type choice. The calculated NBLs are reported by redox facies, if this option has been chosen and at GWB scale. The calculation path type used (A, B, C, D, Table 1) and the confidence level of the resulting NBL (Table 2) are also specified. The results of the calculation as well as the settings used are listed in the second page of the procedure (Figure 2). In addition, the result page provides the links to the intermediate products of the calculation such as the validated dataset (built after time filtering, non-detected values handling and electrical balance verification), pre-selected and facies separated datasets, final datasets used for the calculation and the path type choice. The calculated NBLs are reported by redox facies, if this option has been chosen, and at GWB scale. The calculation path type used (A, B, C, D, Table 1) and the confidence level of the resulting NBL (Table 2) are also specified.

Study Area
Apulia is a southeastern region of Italy ( Figure 3) that covers an area of about 19.500 km 2 hosting more than four millions inhabitants [31]. This region faces the Adriatic and Ionian seas with a coastline extension of about 900 km.

Study Area
Apulia is a southeastern region of Italy ( Figure 3) that covers an area of about 19.500 km 2 hosting more than four millions inhabitants [31]. This region faces the Adriatic and Ionian seas with a coastline extension of about 900 km.
Water 2021, 13, x FOR PEER REVIEW 7 of 21 Figure 3. Map of the Apulia region, showing the main geographical sub-regional areas (Gargano, Tavoliere, Murgia and Salento) and the regional groundwater monitoring network.
Apulia region is characterized by dry sub-humid (Mediterranean) climate with hot and dry summers and mild, wet winters [32]. Precipitation are on annual average around 600 mm, and mainly concentrated between November and February.
The region economy is mainly driven by the agricultural sector followed by service and industry. Agricultural areas cover more than 70% of the region and the main crops are wheat, olive, grapes, citrus, and vegetables.
Apulia region mostly corresponds to the exposed area of the Apulian carbonate platform (Apulia foreland), separated from the Southern Apennines by the Bradanic Through (foredeep). The main geological complexes are the Mesozoic carbonate platform and the Paleocene to Pleistocene sedimentary covers [33] (Figure 4). showing the main geographical sub-regional areas (Gargano, Tavoliere, Murgia and Salento) and the regional groundwater monitoring network. Apulia region is characterized by dry sub-humid (Mediterranean) climate with hot and dry summers and mild, wet winters [32]. Precipitation are on annual average around 600 mm, and mainly concentrated between November and February.
The region economy is mainly driven by the agricultural sector followed by service and industry. Agricultural areas cover more than 70% of the region and the main crops are wheat, olive, grapes, citrus, and vegetables.
Apulia region mostly corresponds to the exposed area of the Apulian carbonate platform (Apulia foreland), separated from the Southern Apennines by the Bradanic Through (foredeep). The main geological complexes are the Mesozoic carbonate platform and the Paleocene to Pleistocene sedimentary covers [33] (Figure 4).  The sedimentary covers consist of heterogeneous marine and alluvial deposits of different ages [34].
All the calcareous and calcareous-dolomite rocks forming the Mesozoic succession are influenced by karst phenomena of different evolution degree. This causes an almost total absence of surface water and a widespread karst landscape. Karstic processes combined to the presence of fractures affecting the Mesozoic succession, give rise to a secondary permeability, which is always markedly discontinuous in space and may vary from very low to high [35,36].
For all of the above, the groundwater resources hosted in the three Mesozoic hydrogeological structures of Gargano, Murgia and Salento, representing the largest coastal karst Italian aquifers, are the main regional water supply of Apulia [37]. They are exposed to contact with the sea and the groundwater systems are largely governed by the freshwater-saltwater equilibria.
The hydrogeological structure of the Gargano is generally characterized by confined groundwater flow, except along a narrow coastal belt. Along the south-western boundary, groundwater interacts with warm waters rising from the deepest part of the carbonate platform [38]. These deep groundwater fluxes are drained by springs in the SW portion of the promontory (Manfredonia Gulf) and in the NW portion toward Lesina and Varano lakes. Particularly, in these areas they assume typical fossil water characteristics [39].
In the Murgia hydrostructure, the presence at different depths of thick dolomitic layers, practically impermeable, causes the carbonate aquifer to be under pressure and with irregular geometry. Only near the coast, and mainly in the southern part of the Murgia, the groundwater flow becomes unconfined. The watershed for groundwater flow approximately overlaps the highest topographic altitudes along the NW-SE direction and allows to identify the Bradanic and the Adriatic sectors.
The Salento is the most permeable carbonate hydrostructure of the Apulia. Groundwater flows mainly under phreatic conditions and the water table heights reach only 3-4 m above sea level. Due to the Salento geographical shape, the groundwater assumes a lenticular shape, floating on seawater. The aquifer mainly is recharged by rainfall, but in the north-west sector, it is also laterally fed by groundwater coming from the Murgia aquifer [36]. The sedimentary covers consist of heterogeneous marine and alluvial deposits of different ages [34].
All the calcareous and calcareous-dolomite rocks forming the Mesozoic succession are influenced by karst phenomena of different evolution degree. This causes an almost total absence of surface water and a widespread karst landscape. Karstic processes combined to the presence of fractures affecting the Mesozoic succession, give rise to a secondary permeability, which is always markedly discontinuous in space and may vary from very low to high [35,36].
For all of the above, the groundwater resources hosted in the three Mesozoic hydrogeological structures of Gargano, Murgia and Salento, representing the largest coastal karst Italian aquifers, are the main regional water supply of Apulia [37]. They are exposed to contact with the sea and the groundwater systems are largely governed by the freshwater-saltwater equilibria.
The hydrogeological structure of the Gargano is generally characterized by confined groundwater flow, except along a narrow coastal belt. Along the south-western boundary, groundwater interacts with warm waters rising from the deepest part of the carbonate platform [38]. These deep groundwater fluxes are drained by springs in the SW portion of the promontory (Manfredonia Gulf) and in the NW portion toward Lesina and Varano lakes. Particularly, in these areas they assume typical fossil water characteristics [39].
In the Murgia hydrostructure, the presence at different depths of thick dolomitic layers, practically impermeable, causes the carbonate aquifer to be under pressure and with irregular geometry. Only near the coast, and mainly in the southern part of the Murgia, the groundwater flow becomes unconfined. The watershed for groundwater flow approximately overlaps the highest topographic altitudes along the NW-SE direction and allows to identify the Bradanic and the Adriatic sectors.
The Salento is the most permeable carbonate hydrostructure of the Apulia. Groundwater flows mainly under phreatic conditions and the water table heights reach only 3-4 m above sea level. Due to the Salento geographical shape, the groundwater assumes a lenticular shape, floating on seawater. The aquifer mainly is recharged by rainfall, but in the north-west sector, it is also laterally fed by groundwater coming from the Murgia aquifer [36].
The fourth main hydrogeological structure of the Apulia Region is the Tavoliere, a shallow aquifer consisting of a complex alternation of different grain size alluvial sediments. This aquifer is unconfined in the western part of the plain, where the coarse-grained sediments prevail, and confined in the central-eastern part where fine-grained material tend to prevail in the upper part of the sequence [40].
Following the criteria established by Italian law for defining the GBWs [5], three main type of hydrogeological complexes (HCs) have been recognized in the Apulia: carbonate, detrital, and alluvial [41,42]. Moreover, on the basis of hydrogeological features of the formation hosting groundwater and its physical-chemical characteristics different GWBs were distinguished within some aquifer for a total of 29 GWBs (Figure 4).
Since the first regional groundwater monitoring network was established in 1995, different redesign have been done until 2015 to fit the WFD provisions and optimize the network configuration [43]. It follows that the surveys have been affected by some discontinuities in time. Currently, the qualitative groundwater monitoring network consists of 341 wells and 20 springs distributed almost uniformly over the regional territory ( Figure 3).
The dataset used in this study consists of more than 3600 measured values of the main hydrogeochemical parameters, discontinuously measured from 1995 to 2018. In general, two sampling campaigns were carried out per year, one at the end of the recharging period and the other at the end of the dry season, for a total of 22.
Since the measurements refer to a very long time interval, the sampling and analytical methods could have likely changed over time. Inhomogeneities are therefore possible, but the NBL assessment procedures have been designed to minimize their effect on the results (use of median values, outliers removal, statistical distribution analysis). Detailed information is available only for data collected after 2000. From this date on, official monitoring protocols and field forms indicate the physical parameters (pH, Eh, EC, temperature and DO) have been measured onsite by a multiparametric probe. Water samples have been filtered through 0.45 µm pore-size filters immediately upon collection and the aliquot for cations and trace elements has been acidified with supra pure HNO 3 1%. All the chemical parameters have been measured and validated at certified public and private laboratories by ion chromatography for anions and by Inductively Coupled Plasma-Mass Spectrometry for major cations and trace elements.
Severe impacts on regional GWBs are due to overexploitation and qualitative degradation [41,42]. Quantitative pressures produce impacts in terms of piezometric lowering, reduction of springs outflow and decrease of groundwater contribution to surface aquatic ecosystems. The impact from quantitative significant pressures can result in an increase of saline intrusion with a consequent reduction of freshwater availability and soil salinization in agriculture.

GW Geochemistry and Descriptive Statistics
A summary of the qualitative status of the Apulia GWBs is provided, based on selected monitoring data.
In Table 3, the main statistics of the physical-chemical parameters relating to the Apulian monitoring network are represented. Groundwater shows a wide range of pH values, with conditions ranging from acidic to frankly alkaline, and slightly alkaline mean values. Oxidizing conditions predominate in the region, although reducing environments are also found within some GWBs, highlighted by DO deficiency and negative ORP values. The variability from 0.1 to tens mS/cm, of the electrical conductivity (EC) is typical of wide aquifers that spread from inland to the coastline.  Table 4 reports the exceedances, with respect to the values prescribed by the Italian water-related Legislative Decrees [4,5], of those parameters resulting as main indicators of significant anthropic pressures and related impacts arisen in the regional GWBs characterization [44,45]. The most evident result is the high percentage of exceedances for chloride (49.0%) and sulphates (16.8%) that can be related to saline contamination phenomena, especially in coastal areas. The percentages relating to nitrate (20.3%) and ammonium (6.4%), classical indicators of anthropogenic contamination, can be linked to the extensive agricultural land use and, partially, to livestock activities. More than a quarter of the dataset shows exceedances of the TVs for iron (32.2%) and manganese (26.2%). Both these parameters are sensitive to the redox reducing conditions of the aquifers which are quite frequent, as mentioned previously. Several values exceeding the legal limits were found for selenium (3.7%) and fluoride (3.5%), while sporadic exceedances were found for boron (0.8%) and arsenic (0.3%).   In more detail, Figure 5 shows the Piper's classification diagram distinguished by HC. This type of representation provides information on the geochemical facies and therefore on the processes that led to the observed main chemistry. Each point represents the median value of Piper's diagram characteristic parameters in each considered GWB, distinguished by color.
Groundwater of the Gargano and Murgia-Salento carbonate HCs (Figure 5a) shows a clear transition from frankly calcium-bicarbonate terms to chloride-alkaline terms, with increasing salinity. The points of the detrital HCs (Figure 5b) show more variable compositions, generally ranging between the chloride-alkaline and the sulphate-chloride-earth alkaline facies. Finally, groundwater from the Miocene carbonate HCs has a clear calciumbicarbonate composition, while the alluvial HCs are characterized by mixed terms, without dominant ions (Figure 5c). The presence of clear bicarbonate-earth alkaline terms is mainly attributable to water-rock interaction processes in the carbonate aquifers, while the enrichment in chlorides, sulfates, and alkalis (sodium in particular), reflecting on a substantial increase in groundwater salinity, are found in coastal GWBs.
The NBLs at the GWB scale were assessed for all the parameters in Table 4, except for nitrogen species, which were used as anthropic indicators. To ensure a certain significance in the statistical processing, the calculation was carried out for the GWBs with at least 7 monitoring points, for a total of 15 investigated GWBs and 294 MSs.

Preliminary Analysis on the Monitoring Stations
The proposed methodology requires some preliminary processing in order to rearrange the dataset for the following steps. At first, all the analytical values below their own limit of quantification (LOQ) must be replaced with a value equal to half the limit itself. Nearly 6500 values have been replaced with regard to the 8 investigated parameters and the 15 considered GWBs.
Successively, the dataset is treated by removing the analyses whose electrical balance error (EBE) exceeds the given threshold. This validation step resulted in the removal of 55 analyses from the overall dataset, whose EBE was exceeding the threshold of 10%.
A third important data preprocessing step concerns the redox facies separation. In the case at hand, once set an ORP threshold of 50 mV (selected following a comparison with the values of Fe and Mn, tendentially showing higher concentrations below this limit), 154 MSs out of 294 were identified in the oxidizing facies and 63 MSs in the reducing one, with distribution and numerical consistency different by GWB. Finally, due to the lack of ORP and DO measurements, the redox facies assignment has not been practicable in 77 MSs. The facies separation was applied only for iron and manganese, which are particularly affected by reducing conditions that increase the solubility of their oxy-hydroxides. Therefore, only for these two metals, the NBLs were first calculated separately for the two redox facies, then a unique value (the highest of the two) was chosen from them. For all the other parameters the NBLs were calculated from the total dataset of 294 MSs.
The last step of the preliminary analysis is the pre-selection of the MSs. In this phase, some MSs and the related datasets are removed based on an anthropogenic contamination test. The test operates using two different markers of anthropogenic contamination, according to the redox conditions: a threshold of 37.5 mg/L of nitrate for the oxidizing facies and 0.375 mg/L of ammonium for the reducing one. Only those MSs in oxidizing/reducing facies, where nitrates/ammonium do not exceed the respective concentration threshold, have been retained for next processing; all the others have been considered somehow contaminated by anthropic activities and discarded. The selected concentration limits correspond to 75% of the expected quality standards/TVs. When the methodology does not involve the facies separation, both markers were used for the MSs pre-selection. On average, 66% of MSs complies with the adopted pre-selection criteria, but with very diversified distributions within the different GWBs.
In particular, the pre-selection criterion allowed retaining a number of monitoring stations (pre-selected monitoring sites = PS MSs) ranging from 40% to 90% of the initial number in 10 of the considered GWBs. The percentage of pre-selected stations falls down to values smaller than 38% in the remaining 5 GWBs. Unfortunately, these latter five GWBs, mostly belonging to the Tavoliere aquifer, were initially affected by a poor number of MSs, so the related computed NBLs are expected to be characterized by a very low confidence level. In conclusion, 194 monitoring sites were selected for the NBLs assessment in 15 GWBs of which 104 in the oxidizing and 54 in the reducing facies (Table 5).

NBLs Calculation for the Non-Redox Sensitive Elements
For each PS MSs, a representative value was defined as the median of all available measurements for each investigated parameter (fluoride, chloride, sulphate, boron, arsenic, selenium). The outliers were then analyzed using the Huber's non-parametric test and the identified anomalous data were eliminated from the dataset, before performing the calculation of the NBLs. Consequently, the number of the representative values is not equal for all parameters, resulting in part from the elimination of the outliers and in part from any missing analysis for that specific parameter.
According to the guidelines, the specific path types (A, B, C, D) were identified based on the pre-selected dataset size. For type A and B the Shapiro-Wilk test was then applied to verify the normality of the dataset. Data relating to the performed elaborations (including the specific path types and the number of representative values for each parameter) are shown in Table S1 (Supplementary Material).
The most frequent type of dataset were C and D, then A and finally B were the least represented (see Table S1 in the Supplementary Materials).
The defined NBLs for all the parameters at the GWBs scale and the associated confidence levels defined considering the aquifer type (confined or unconfined), its extension and the number of data used for the NBL assessment, are shown in Table 6 and in Figure 6. Chlorides show diffusely high background values, often above the corresponding TV (80% of the selected GWBs). Values sometimes higher than the legal limits have also been calculated for sulphate (6 GWBs) and selenium (4 GWBs), while the NBLs defined for fluoride, boron and arsenic are always lower than the respective TVs. Unfortunately, only for 4 GWBs it was possible to associate a high confidence level to the calculated NBLs, while it resulted almost always low or very low. The low/very low confidence levels are attributable to a strong reduction of MS number due to the high percentage of MS exceeding nitrate/ammonia and/or the lack of the specific analytical value.

NBLs Calculation for the Redox Sensitive Elements (Iron and Manganese)
As already mentioned, for iron and manganese the NBLs were calculated separately for the oxidizing and reducing facies, and then the higher of the two values was chosen as NBL of the entire GWB. When the number of data did not allow the NBL calculation for both facies, the only available value was assigned to the GWB. The calculated NBLs and the associated confidence levels are shown in Table 7 and Figure 7, while the other data relating to the performed elaborations are shown in Table S2 (Supplementary Materials).

NBLs Calculation for the Redox Sensitive Elements (Iron and Manganese)
As already mentioned, for iron and manganese the NBLs were calculated separately for the oxidizing and reducing facies, and then the higher of the two values was chosen as NBL of the entire GWB. When the number of data did not allow the NBL calculation for both facies, the only available value was assigned to the GWB. The calculated NBLs and the associated confidence levels are shown in Table 7 and Figure 7, while the other data relating to the performed elaborations are shown in Table S2 (Supplementary Materials).   For both elements, the exceedances of the TVs concern almost all GWBs and are mainly due, as expected, to values defined for the reducing facies. Particularly, iron shows severe exceedances, with NBLs even two orders of magnitude higher than the corresponding TV (200 µg/L).
The confidence levels associated with the defined NBLs are always low or very low. This is due to the reduced number of initial MSs (217), as well as to the redox facies separation, which produces two different datasets with reduced numerical consistency.

Discussion and Conclusions
The Apulian GWBs characterization, carried out according to the rules set by the national regulation [5], pointed out the exceedance of the threshold values for some inorganic parameters, identifying a group of groundwater bodies at risk of failing to meet the Water Framework Directive environmental objectives. Consequently, the assessment of the natural background levels becomes an urgent need, in order to establish the natural or anthropogenic origin of the high values encountered for such parameters.
This work reports a tentative application of the recently published national guidelines for the NBL assessment to the inorganic compounds monitored in the GWBs of the Apulia region. The NBLs values calculated resulted often higher than the related threshold values.
Concerning the redox-sensitive elements iron and manganese, the estimated NBLs exceed, even of 11 and 10 orders of magnitude respectively, the TVs at 15 GWBs. The lack of redox-related parameter measurements (i.e. ORP and DO) forced us to eliminate many MSs from the available database, making the pre-selected dataset rather small and the estimated NBLs confidence levels low or very low. Furthermore, the resulting NBLs almost always correspond to the values estimated in the reducing facies confirming the extreme sensitivity of these two elements in a negative redox environment. This outcome

Discussion and Conclusions
The Apulian GWBs characterization, carried out according to the rules set by the national regulation [5], pointed out the exceedance of the threshold values for some inorganic parameters, identifying a group of groundwater bodies at risk of failing to meet the Water Framework Directive environmental objectives. Consequently, the assessment of the natural background levels becomes an urgent need, in order to establish the natural or anthropogenic origin of the high values encountered for such parameters.
This work reports a tentative application of the recently published national guidelines for the NBL assessment to the inorganic compounds monitored in the GWBs of the Apulia region. The NBLs values calculated resulted often higher than the related threshold values.
Concerning the redox-sensitive elements iron and manganese, the estimated NBLs exceed, even of 11 and 10 orders of magnitude respectively, the TVs at 15 GWBs. The lack of redox-related parameter measurements (i.e., ORP and DO) forced us to eliminate many MSs from the available database, making the pre-selected dataset rather small and the estimated NBLs confidence levels low or very low. Furthermore, the resulting NBLs almost always correspond to the values estimated in the reducing facies confirming the extreme sensitivity of these two elements in a negative redox environment. This outcome is not surprising considering that Apulian GWBs, and particularly the carbonate ones, are often very deep and confined with the roof at even more than 100 m under the sea level.
The abundance of iron and manganese in the carbonate Apulian GWBs can be associated both to clay-rich lenses interbedded in the carbonate sequence and to the widespread red soils covering the limestone and filling karst cavities and fractures. Their presence is also justified in some detrital and alluvial aquifers that derive from the weathering and erosion of carbonate rocks [44,45].
The estimated NBLs of chloride and sulphate exceed the corresponding TVs in 12 and in 6 of the 15 investigated GWBs respectively, with a wide range of variability. This result was easily predictable for the coastal GWBs but unexpected in the Murgia Bradanica, being this latter very far from the coastline. As known, chloride concentration in coastal GWBs can be particularly affected by the current seawater intrusion, nevertheless different concurrent saline sources can overlap to seawater contributing to groundwater salinization independently from the GWB distance from the coastline [26,27].
In some cases, the mobilization of paleo-seawaters intruded into the aquifers during previous marine transgressions can occur. These saline waters, isolated from active flow, reside for a period long enough to allow water to assume typical connate water characteristics. This is the case of the deep groundwater flux exchanges that occur along the main tectonic elements as an effect of the Apennine convergence towards the Apulia foreland giving rise to upward flow of warm waters drained by springs [38,39]. Previous studies revealed that the chemical and isotopic features of salt groundwater, sampled in different point of the Apulian karstic aquifers or drained from coastal springs, are very different from modern seawater [46][47][48][49]. Nevertheless, currently it is not clear to what extent this phenomenon can be considered natural or forced by anthropogenic activity (e.g., excessive pumping in coastal area).
In 4 of the 15 considered GWBs, the calculated NBL for selenium exceed the threshold value set for this element. Selenium can be present in trace in the carbonate rocks and the argillaceous sediments that fill up the karst cavities [50,51]. In many cases, high values of selenium concentration have been measured in springs draining deep water whose composition results from a high degree of water-rock interaction. Nevertheless, given the evident spatial and temporal fragmentation of high selenium concentration measurements, some Authors proposed the hypothesis that it may originate from its use in agriculture as a nutrient, which, in particular environmental circumstances (e.g., soil temperature, humidity and acidity), leaches from the ground into the underlying groundwater [52].
The application of the Italian Guideline to the estimation of the NBLs of the Apulian GWBs, as proposed in this work, highlighted some problems related to the procedure itself. The Guidelines require an initial preselection of the monitoring sites and the related sets of analysis to retain for the NBLs estimation, based on several validation tests. Unfortunately, this preliminary step often reduces consistently the available information making the estimations unreliable or impossible. This is particularly true in those groundwater bodies characterized by a bounded extension and usually monitored with a few sampling sites. It is not a case that the NBLs estimation has not been possible in almost half of the regional groundwater bodies.
The data spatial consistency and the temporal continuity of monitoring are critical, but the existing groundwater monitoring networks have been designed for different purposes, such as the evaluation of the environmental state of the groundwater body or keeping under control local-specific issues.
The used dataset covers a time span of 24 years during which the monitoring has been suspended even for years and then restarted, sampling and analytical methods have changed with the scientific and technical advance. As a consequence, spatial consistency and temporal continuity of the available data can result deficient concerning some parameter or the whole groundwater body.
All the above requires a deep technical and economical evaluation of the opportunity of redesigning or modifying the existing monitoring network in order to reduce the weakness highlighted with regard to the natural background levels estimation in the regional groundwater bodies. This topic is the subject of specific further investigations in an experimental coastal area, whose final goal is to solve the key question of source of salinization, and define specific protocols for NBLs assessment in coastal aquifers [53][54][55].
In conclusion, the aim of this work was to apply and verify the national guidelines for the assessment of the NBLs, which implement the EU Water Framework Directive, in a south-European, regional contest characterized by the typical Mediterranean climatic and hydrologic features. A frequent issue of such areas is the high anthropization of the coastal areas and the overexploitation of groundwater for agricultural needs in almost total lacking surface resources. Even though some local specific issues arose during the application of the methodology, the results are suitable to be applied for the next steps of the EU-WFD implementation. In general, the NBL definition can be considered, at any latitude, one of the pillars for sustainable and long-term groundwater management, by tracing a clear boundary between what is supposed to be natural and what is clearly anthropic.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/w13070958/s1, Table S1: Results of the procedure for the NBLs calculation at the GWB scale, applied to the non-redox sensitive elements, Table S2: Results of the procedure for the NBLs calculation at the GWB scale, applied to the redox sensitive elements (Fe and Mn). Data Availability Statement: Paper experimental data are not currently publicly available but can be supplied on request from the corresponding author.