Numerical Groundwater Model to Assess the Fate of Nitrates in the Coastal Aquifer of Arborea (Sardinia, Italy)

Abstract

The Arborea plain in Sardinia (Italy) is classified as a nitrate vulnerable zone (NVZ). In the present study, the individual work steps that are necessary to progress from the existing 3D hydrogeological model to a 3D numerical groundwater model using the interactive finite-element simulation system FEFLOW 7.4 are shown. The results of the transient flow model highlight the influence of the drainage network on the overall groundwater management: the total water volume drained by the ditches accounted for approximately 58% of the annual outflow volume. The numerical transport simulations conducted from 2012 to 2020 using hypothetical field-based nitrate input scenarios globally underestimated the high concentrations that were observed in the NVZ. However, as observed in the field, the computed nitrate concentrations in December 2020 still varied strongly in space, from several mg L⁻¹ to several hundreds of mg L⁻¹. The origin of these remaining local hotspots is not yet known. The modeling of rainfall fluctuations under the influence of climate change revealed a general long-term decline in the groundwater level of several tens of centimeters in the long term and, in conjunction with a zero-nitrate scenario, led to a significant decrease in nitrate pollution. Although hotspots were attenuated, the concentrations at several monitoring wells still exceeded the limit value of 50 mg L⁻¹.

1. Introduction

Water resource management and water governance in Mediterranean (MED) coastal zones are important and urgent challenges. The development of large cities and megacities is putting increasing pressure on the quality and quantity of water resources because of the increasing population and the concentration of economic activities in MED coastal zones [1,2]. Coastal areas around the Mediterranean basin are often highly populated and concentrated, with multisector economic and agricultural activities. This has resulted in an important need for fresh water and a high solicitation of coastal aquifers, which can lead to saltwater intrusion [3]. This issue, combined with contaminated surface water that percolates towards these aquifers, as well as climate change, reveals the need for innovative groundwater management, especially in coastal areas [4].

Numerical modeling is one of the most relevant tools that can be used to investigate and understand groundwater systems, estimate the effects of certain groundwater management systems, reconstruct the evolution of a hydrological system over time, and forecast the effects of future actions or hydrological and climatic evolutions. In order to achieve these results with numerical modeling, it is essential to keep in mind that a model, as it may be articulated, cannot represent the real world in every detail; a model is always a simplified representation of the complex reality of a natural system [5]. The mathematical modeling of groundwater flow and water balances in alluvial aquifers, therefore, requires time-dependent input data such as water fluxes within the modelled region and across its boundaries [6,7,8,9,10]. While in many field studies, local pumping and infiltration rates are most accurately determined as an element of the global water balance, the quantification of other components, such as natural groundwater recharge from rainfall and lateral subsurface inflow, may be rather difficult.

Groundwater recharge due to precipitation is an important part of any water balance evaluation. Although it is one of the most relevant components in groundwater studies, recharge is also one of the least understood since its rates can vary widely over time and in space and are difficult to measure directly [11]. Recharge depends on numerous factors, such as precipitation, evaporation, transpiration, surface runoff, and climatic conditions at the soil surface, making it difficult to accurately quantify the infiltrated water flux at the water table [11,12,13,14,15,16,17,18]. Furthermore, the soil lithology, soil texture, and, particularly, the depth of the water table play critical roles in groundwater recharge [19,20,21,22,23].

There are no direct methods for quantifying subsurface inflow and outflow. A descriptive evaluation of the subsurface inflow rate can be made using environmental isotopes [24,25,26,27]. To obtain an estimation of the subsurface flow, a water flux per meter boundary is usually calculated by applying Darcy’s law based on both the measured hydraulic gradient across the boundary and the transmissivity of the aquifer. However, quite often in numerical studies, the flow rates at these boundaries, commonly called prescribed flux boundaries, are obtained by calibration of the groundwater flow model. Sometimes, model users apply a water balance of the catchment area near the chosen boundary to obtain estimates of boundary fluxes [28,29,30] and, depending on the model approach used, an approximation of the delay of lateral inflow due to precipitation events.

Four case studies located at the two shores of the MED Sea (Greece, Turkey, Tunisia, and Italy) are part of the PRIMA project Sustain-COAST. In the present paper, one of these study areas, the coastal aquifer system of the Arborea plain in Sardinia (Italy), is analyzed. It is characterized by intense agricultural activity based on dairy cattle farming. The area was reclaimed from a lagoon approximately 100 years ago. Thus, an important drainage network has been developed to maintain the soil in a suitable condition for agriculture. Its aquifer system is highly contaminated by nitrates. The Arborea plain has been classified as a nitrate vulnerable zone (NVZ) since 2005 [31,32].

Although nitrogen (N) is essential for healthy plant and animal populations, high concentrations of this nutrient can degrade water quality through eutrophication. Excessive concentrations of nitrate (NO₃) (the most common form of N dissolved in streams and groundwater) in surface water can stimulate the growth of algae and nuisance organisms. A significant part of the total N load in streams originates from groundwater discharge to surface water, open-air ditches, and the sea. This is particularly true for the Arborea plain, where the base flow is a dominant component of the dense drainage network of open-air ditches. Previous research based on the combination of multi-scale hydrogeological and hydrochemical field studies, groundwater flow and transport modeling, dual stable isotope approaches, and geochemical modeling has mainly indicated agriculture and inappropriate sewage water management as the sources of NO₃ contamination of the groundwater, which, moreover, is also affected by geochemical processes [33].

To reduce and prevent the water pollution caused or induced by NO₃ from agricultural sources, the Nitrates Directive [34] requires the Member States to designate NVZs, i.e., areas where surface waters contain an excessive concentration of NO₃, that is, where the groundwater NO₃ levels are likely to or already exceed 50 mg L⁻¹. Member states were required to limit the amount of manure N that is applied to the land to 210 kg N ha⁻¹ yr⁻¹ by 1999 for each farm or livestock unit. In 2003, for example, this value was reduced to 170 kg N ha⁻¹ yr⁻¹.

The Arborea plain represents a relevant case study in Sardinia, where several research activities have been performed in the last 15 years through an interdisciplinary approach that includes geology, hydrogeology, hydrogeochemistry, and agronomy to assess the process of groundwater contamination by nitrate. An agronomic study at the field scale was carried out to evaluate the impact of four fertilization systems on nitrate losses in a double-crop system (maize and Italian ryegrass) [35]. Through monthly monitoring of the nitrate concentrations in the soil solution (2009–2012), clear seasonal dynamics of the nitrate concentration were observed, with the maximum occurring in autumn–winter. The main findings of this study showed the following: (i) replacing organic sources of N with mineral ones did not substantially reduce nitrate leaching; (ii) nitrate leaching was unavoidable in autumn–winter; (iii) the nitrate concentrations were still high in late summer despite the high N removal of the maize crop; and (iv) the minimum release of nitrate during the leaching period was observed when using only a high organic C/N ratio as fertilizer for the ryegrass, but yield reductions were also observed, especially for the ryegrass. A hydrogeological and hydrogeochemical study was carried out over the scale of the whole plain to reconstruct a 3D hydrogeological model and define the major groundwater flow paths and the distribution of nitrate concentrations in the aquifers [36]. Several highly polluted areas were identified in the northeast sector of the NVZ and in the southern part of the plain, where an increasing trend from east to west in the direction of the main groundwater flow prevailed; these results indicated an additional pollution source located outside the NVZ. Hydrogeochemical and multi-isotopic studies were performed by combining environmental isotopes (δ¹⁵N–NO₃, δ¹⁸O–NO₃, δ³⁴S–SO₄, δ¹⁸O–SO₄, δ¹³C–DIC, δ¹¹B), water quality, and hydrogeological indicators to assess the nitrate sources and their fate in the groundwater as well as the occurrence of nitrate attenuation processes such as denitrification. The nitrate isotopic composition confirmed the occurrence of denitrification processes, which were consistent with other geochemical indicators, such as low dissolved oxygen concentrations and low Eh values, indicating conditions suitable for denitrification in the aquifer [37]. The integration of the δ¹¹B data with δ¹⁵N values allowed for the assessment of the nitrate source in the study area, clearly showing that organic fertilizers were the main source of nitrates in the groundwater [32].

Some modeling applications were carried out in the Arborea area. Cau and Paniconi [38] applied a semi-distributed Soil and Water Assessment Tool (SWAT) to predict the impact of different land management practices on water and the agricultural chemical yield over a long period. An analysis matrix was developed by combining water quality and quantity indicators for the rivers, lagoons, and soil with socioeconomic variables using a multicriteria decision support system. The results suggested that the most widely acceptable option consisted of the transfer of intensive agricultural practices to the larger watershed, which is less vulnerable, in tandem with the reuse of treated wastewater for irrigation purposes. Matzeu [39] estimated the nitrate concentrations in groundwater by using numerical methods; a numerical transport model was obtained by applying the 3D Multi-Species Transport Model (MT3D) and coupling it with a groundwater flow model developed by using the Modular 3D Finite-Difference Groundwater Flow (MODFLOW) model. The results showed that higher nitrate values were located in areas characterized by the highest agricultural and livestock farming concentrations.

The need for the improvement of the existing groundwater flow model of Matzeu [39] was identified mainly through three physical aspects. These concerned the following: (i) the implementation of a more hydrologically oriented representation of the inflow rates resulting from adjacent catchments at the eastern model boundary, (ii) the evaluation of spatially distributed unsteady groundwater recharge in a simplified manner, and (iii) the consideration of the given dense drainage network of open ditches. For example, Matzeu [39] used the general head boundary condition to evaluate the unsteady inflow rate at the eastern boundary. The inflow rate was the difference between an assumed constant hydraulic head along the boundary and the groundwater head simulated by the numerical model multiplied by an estimated value of the conductance along the boundary. Regarding groundwater recharge, Matzeu [39] derived the transient variation in the vertical inflow rate in the groundwater from a water balance analysis of the unsaturated zone.

The methodology presented here is based on a novel approach to modeling both the transient groundwater flow in the coastal area of the Arborea plain, taking into account hydrology-based estimated lateral flow rates on the eastern boundary, zonal groundwater recharge using potential infiltration coefficients depending on the soil texture at the soil surface, and the main features of the drainage network, and the fate of leached NO₃ from the unsaturated zone into the aquifer system using the numerical tool finite-element subsurface FLOW (FEFLOW) 7.4. The present study implemented improved datasets for hydrological modeling, a physiography-based water balance to estimate the aquifer recharge that may be applied to other groundwater systems, and field-based nitrate input scenarios to assess the actual state of pollution.

The objective of this study was to validate the developed numerical flow and transport model for the prediction of the water quantity and water quality in the study area. Different scenarios of nitrate leaching were compared to assess the actual fate of the nitrate pollution in the sandy hydrogeological unit of the given complex aquifer system. The modeling approach included climate change scenarios to model the influence of future rainfall on the predicted groundwater levels, which may have a significant impact on both available groundwater reserves and groundwater quality.

In this regard, this study aimed to address the following questions:

• Does the given drainage network have a significant impact on the water balance? If so, to what extent?
• How can the lateral inflow rates be accurately implemented on a hydrologically based assumption?
• The application of nitrate directive prescriptions to agricultural practice in the Arborea plain have not provided satisfactory results so far, as shown by the trimestral monitoring carried out in wells and piezometers by the local environmental agency. Can 9-year simulations made with a numerical model based on relatively low-input and calibration datasets that take into account field-based nitrate input scenarios help in diagnostics to design and implement effective mitigation strategies?
• Does the developed model visibly predict a significant reduction in nitrate mass in the groundwater, or are the observed highly polluted areas also predicted by the model?

2. Materials and Methods

2.1. Study Area Description

The study area is located in the northern part of the Campidano plain (central–western Sardinia, Italy). The volcanic complexes of Montiferru and Monte Arci represent the northern and eastern borders, respectively. The area is limited to the south by the Rio Mogoro River and the Marceddì and San Giovanni Lagoons and to the west by the Oristano Gulf. Approximately 60 km² of the plain is occupied by the Arborea farming district, which is located between the coast and the reclaimed Sassu Lagoon (Figure 1).

The climate is Mediterranean, with a mean annual temperature and precipitation (concentrated between October and March) of 16.7 °C and 575 mm, respectively, and an annual ET₀ of 1164 mm, corresponding to an aridity index of 0.49 [40].

2.1.1. History

The study area was an insalubrious swamp infested by mosquitoes that spread malaria. In the 1920s and 1930s, extensive land reclamation work was implemented for the entire plain: sand dunes were flattened, and brackish and salted wetlands were drained by pumping water from below sea level. Fertile agricultural land was generated, organized in rectangular fields of 2 to 4 ha in size and surrounded by a drainage network consisting of main channels and a dense network of smaller drainage channels with E-W and N-S directions. This drainage system also includes some dewatering pumping stations for water flow regulation in the main channels, such as those of Sassu and Luri, which pump approximately 75 hm³ yr⁻¹ and 4 hm³ yr⁻¹ of drainage water on average, respectively. The reclaimed land was colonized in the 1930s by peasants from northeastern Italy, who formed a very strong cooperative system. Currently, the farmers’ cooperative includes more than 200 farms that manage approximately 30,000 dairy cattle on a 6000 ha irrigated plain, representing one of the most productive agricultural sites in Italy [35].

This intensive dairy cattle system has caused nitrate groundwater pollution that originated from the intensive input of effluents (slurry and manure). Consequently, the Arborea area was identified as an NVZ in 2005 [36]. The forage cropping system for dairy livestock is currently based on the double-crop rotation of silage maize (Zea mays) and Italian ryegrass (Lolium multiflorum), representing over 80% of the irrigated plain. Almost all the remaining area in the district is covered with meadows of alfalfa (Medicago sativa) and horticultural crops (watermelon, melon, strawberries, etc.).

2.1.2. Hydrogeology

As reported by Ghiglieri et al. [36], three hydrogeological units have been identified in the Arborea plain: a sandy hydrogeological unit (SHU), an alluvial hydrogeological unit (AHU), and a volcanic hydrogeological unit (VHU) (Figure 2). The SHU hosts a phreatic aquifer in the Holocene littoral sands cropping out in the plain. The aquifer is characterized by good permeability, with a K value ranging between 10⁻⁵ and 10⁻⁶ m s⁻¹. The AHU is a multilayer aquifer with good permeability (K = 10⁻⁴–10⁻⁵ m s⁻¹) hosted in Pleistocene continental deposits. This aquifer is confined in the plain because it is separated from the upper SHU by an aquitard of silt and clay that outcrops in the Sassu Lagoon. As this aquitard is lacking in the southern part of the plain, the SHU and AHU aquifers are in hydraulic communication with each other (Figure 2). Finally, the VHU, which is located in the volcanic formations of Monte Arci, was not considered in the modeling because of its characteristics and poor exploitation by few wells along the mountain.

2.1.3. Water Use

The Arborea farming district receives high-pressure irrigation water from the Eleonora d’Arborea dam (some 35 hm³ yr⁻¹), located approximately 35 km away. The Oristanese Land Reclamation Consortium distributes an annual amount of surface water equal to 28.4 hm³ yr⁻¹ at the field scale by a permanent sprinkler system for a total irrigated area of 4900 ha (approximately 5800 m³ ha⁻¹). Groundwater, mainly from the SHU aquifer, is used by farmers only for cattle watering (approximately 0.97 hm³ yr⁻¹) and for the washing of livestock facilities (approximately 0.04 hm³ yr⁻¹) (source: LAORE—Regional Agency of Agricultural Policies Application and Rural Development). The Farmers’ Cooperative and the Milk Processing Cooperative 3A use groundwater from the AHU aquifer for industrial processes in their facilities (0.6 hm³ yr⁻¹ and 0.27 hm³ yr⁻¹, respectively).

2.1.4. Groundwater Quality

Since 2007, ARPAS has monitored the groundwater quality in the NVZ with a tri-monthly frequency. The monitoring network consists of 43 sampling points (27 piezometers and 16 wells). In

Table 1. Average annual concentrations of nitrate, ammonium, and nitrite from the groundwater monitoring network (2007–2020) of ARPAS. Values in bold show average annual concentrations of nitrate, ammonium, and nitrite above the thresholds of 50 mg L⁻¹, 0.5 mg L⁻¹, and 0.5 mg L⁻¹, respectively.

ID_WELL (ARPAS Network)	X Coordinate (WGS 84)	Y Coordinate (WGS 84)	NO3		NH4		NO2
			Average	SD	Average	SD	Average	SD
			mg L−1	mg L−1	mg L−1	mg L−1	mg L−1	mg L−1
17PZ008ar	461,524	4,407,087	11.55	11.78	0.75	0.64	0.02	0.00
P01ar	461,793	4,404,058	12.88	7.80	0.08	0.06	0.43	0.53
P02ar	462,647	4,404,165	171.51	26.79	0.12	0.12	1.38	0.85
P03ar	463,755	4,404,111	23.76	26.02	0.26	0.24	0.30	0.83
P04ar	461,798	4,402,316	329.81	59.81	0.17	0.25	3.06	1.44
P05ar	461,060	4,401,284	32.93	38.13	0.12	0.13	0.63	0.48
P06ar	462,012	4,401,286	98.45	54.44	0.34	0.43	0.72	0.45
P07ar	462,652	4,401,289	30.05	15.32	0.08	0.05	0.32	0.21
P08ar	463,609	4,401,302	97.76	37.43	0.05	0.03	0.18	0.10
P09ar	464,530	4,400,935	83.07	69.89	1.02	1.33	0.24	0.17
P10ar	461,269	4,399,687	125.11	23.07	0.16	0.23	1.94	0.78
P11ar	463,166	4,398,687	115.62	37.07	0.13	0.11	0.97	0.74
P12ar	462,203	4,397,650	189.43	74.92	0.06	0.10	0.09	0.13
P13ar	462,124	4,401,291	76.25	29.89	0.12	0.20	0.90	0.46
P14ar	462,975	4,400,705	181.03	44.61	0.04	0.05	0.56	0.42
P15ar	464,029	4,401,664	89.02	30.46	0.06	0.06	0.31	0.37
P16ar	464,276	4,401,706	43.13	23.78	0.06	0.17	0.42	0.27
P17ar	462,650	4,398,853	44.09	15.71	0.06	0.07	0.78	0.60
P20ar	463,601	4,406,282	286.70	70.92	0.07	0.06	0.60	0.62
P21ar	459,789	4,398,871	22.01	31.33	0.32	0.20	0.72	0.75
P22ar	461,451	4,399,488	148.99	37.82	0.14	0.16	2.07	0.84
P23ar	458,950	4,397,389	198.40	84.22	0.12	0.15	2.13	1.30
P24ar	461,293	4,396,566	57.34	32.41	0.77	1.41	0.48	0.30
P26ar	464,449	4,406,075	226.38	26.00	0.03	0.03	0.04	0.02
P27ar	464,024	4,406,258	64.33	12.57	0.05	0.06	0.42	0.17
P28ar	463,399	4,401,307	80.38	23.61	0.06	0.07	0.32	0.08
P29ar	462,651	4,400,515	2.96	5.22	0.21	0.07	0.12	0.34
P30ar	460,104	4,398,210	3.06	7.81	1.18	0.19	0.15	0.25
P31ar	462,498	4,406,741	2.50	2.75	0.21	0.08	0.04	0.03
P32ar	462,838	4,406,745	1.50	1.09	0.15	0.11	0.03	0.02
P34ar	461,572	4,400,003	53.38	32.38	0.07	0.07	0.17	0.17
P35ar	462,881	4,402,472	68.51	16.49	0.08	0.11	0.98	0.71
P36ar	460,342	4,400,943	40.04	18.60	0.08	0.06	0.43	0.27
P37ar	460,418	4,400,804	5.73	19.21	0.19	0.07	0.09	0.09
P38ar	460,691	4,396,128	141.24	46.03	0.07	0.13	0.12	0.06
P39ar	465,035	4,404,262	56.37	38.81	0.04	0.02	0.05	0.03
P40ar	464,805	4,402,660	53.73	16.24	0.03	0.01	0.10	0.05
P41ar	464,823	4,399,386	32.35	18.71	0.10	0.14	0.13	0.27
P42ar	461,452	4,396,954	56.97	20.58	0.03	0.01	0.10	0.08
P43ar	459,512	4,396,949	84.23	24.77	0.38	0.66	0.36	0.24
P45ar	461,841	4,396,172	3.25	2.58	0.29	0.61	0.10	0.14
P46ar	463,935	4,397,919	5.54	1.91	0.03	0.01	0.06	0.04
P47ar	461,109	4,397,773	14.02	23.53	0.15	0.09	0.25	0.34

, the average annual nitrate, ammonium, and nitrite concentrations of the 43 groundwater sampling points during the monitored period (2007–2020) are reported. As indicated by the values in bold, 25 points out of 43 (58%), distributed in the whole NVZ, showed average annual nitrate concentrations up to 6 times above the threshold value of 50 mg L⁻¹. Regarding the other N compounds, only 4 sampling points (9%) showed average annual ammonium concentrations above the threshold value of 0.5 mg L⁻¹, while 14 of them (32%) had average annual nitrite concentrations up to 6 times higher than the threshold limit (0.5 mg L⁻¹).

2.2. Numerical Groundwater Model

The starting point of this study was the hydrogeological conceptual model developed by Ghiglieri et al. [36]. A 3D numerical groundwater model was built using the interactive finite-element simulation system FEFLOW 7.4 [41]. It is one of the most comprehensive, well-tested, and reliable programs for the simulation of flow, groundwater age, and mass-transport processes in porous media.

The model domain of the study area is characterized by two main aquifers; the upper aquifer, which is unconfined and hosted in the SHU unit, is separated from the second aquifer, which is hosted in the AHU unit, by aquitard lagoon deposits. The aquitard has been shown to become thinner towards the south of the area until both aquifers come into contact [36]. To consider the aquitard discontinuity, we used a partial unstructured mesh. The model used allowed for the prediction of the impact of various management scenarios and anthropogenic activities on the progression of the water quantity and quality in time and space.

2.2.1. Governing Equations of the Numerical Flow and Transport Model

FEFLOW 7.4 is a modular 3D finite-element groundwater flow model capable of simulating transient processes of flow and mass transport in subsurface water resources [41].

In the present study, the solute was considered to be nonreactive, i.e., it was considered to behave like a conservative tracer in the groundwater flow. The basis of modeling the flow and transport of a nonreactive solutes in saturated porous media is given by the fundamental physical principles of the conservation of both momentum and mass at the scale of a Representative Elementary Volume (REV) [42].

The general flow equation for saturated groundwater flow is formulated by applying the law of conservation of mass over a macroscopic control volume of porous media located in the flow field [42]. The net inflow into the volume must equal the rate at which water accumulates within the volume under investigation [9], which leads to

$$S_s{\displaystyle \frac{\partial h}{\partial t}}={\displaystyle \frac{\partial }{\partial x_i}}\left(K_{ij}{\displaystyle \frac{\partial h}{\partial x_j}}\right)-Q$$

(1)

where subscripts i,j (=1,2,3) represent the principal coordinate directions, x [L] represents the space coordinates, t [T] represents the time, K_ij [L T⁻¹] represents the tensor of hydraulic conductivity, h [L] represents the hydraulic head, S_s [L⁻¹] represents the specific storage, and Q [T⁻¹] represents local sources and sinks per unit volume.

The general three-dimensional equation describing the contaminant transport of a nonreactive solute in groundwater is described by the following second-order partial differential equation [9]

$$\begin{array}{c}{\displaystyle \frac{\partial C}{\partial t}}=-{\displaystyle \frac{\partial }{\partial x_i}}\left(v_{i}C\right)+{\displaystyle \frac{\partial }{\partial x_i}}\left(D_{ij}{\displaystyle \frac{\partial C}{\partial x_j}}\right)+D_m{\displaystyle \frac{{\partial }^{2}C}{\partial x_i^2}}-Q{C}_{in}\\ \mathrm{with} v_i=-{{\displaystyle \frac{1}{n_c}}K}_{ij}{\displaystyle \frac{\partial h}{\partial x_j}} ,\end{array}$$

(2)

where C [M L⁻³] represents the concentration of the solute in the water, v_i [L T⁻¹] represents the mean flow velocity, D_ij [L² T⁻¹] represents the dispersion tensor [8], D_m [L² T⁻¹] represents the molecular diffusion coefficient, n_c [-] represents the kinematic porosity, and C_in [M L⁻³] represents the concentration of the sink/source term.

The set of differential equations that allow for the study of contaminant transport in a groundwater system is therefore obtained in the form of a coupled system of two types of equations: a groundwater flow equation (Equation (1)) and a contaminant transport equation (Equation (2)).

The Galerkin finite-element method was adopted to solve these two equations using an automatic time step control via predictor–corrector schemes in the FEFLOW subroutines to perform simulations. Initial and boundary conditions are required for the solution of the flow and transport equations. The initial conditions are the hydraulic head and the concentrations at the starting time of the simulation. Boundary conditions related to the two primary variables (h, C) must be specified for the entire model boundary, and all may vary with time.

2.2.2. Geometry of the Model Domain

The geometry of the numerical model of the Arborea coastal plain was reconstructed on the basis of previous studies in the area and was based on the 3D hydrogeological conceptual model developed by Ghiglieri et al. [36] (Figure 3). The surfaces that delimit the two main hydrogeological units, i.e., the SHU and AHU, and the surfaces that delimit the aquitard were extrapolated from the 3D model of Ghiglieri et al. [36] and were implemented in the numerical model (Figure 4). The digital terrain model (DTM) with a 10 m resolution was provided by the “Sardegna Geoportale” resource of the Sardinia region [43]. This was used as the first surface (or slice n°1) of the partially unstructured model and was the surface at which the outcropping layers, i.e., the SHU, the AHU, and the aquitard, were pinched using the specific tool in the 3D layer configuration dialogue in FEFLOW 7.4.

The model domain was spatially discretized into triangular grids containing 58,022 nodes and 91,798 elements, as shown in Figure 4. Grid sizes were locally refined to set the nodes on the observation wells and to represent the given physical model boundaries. To reduce the numerical effort, only the main branches of the drainage network were considered; the locations of the ditches were approximated by using the closest grid nodes.

The hydrodynamic parameters used for the different aquifer units were chosen in accordance with the hydrogeological model of Ghiglieri et al. [36]. The hydraulic isotropic conductivities of the AHU unit and aquitard were chosen to be constant at 6.048 m d⁻¹ and 0.000864 m d⁻¹, respectively. A spatial variation in the permeability field within the SHU was considered; the optimized distribution obtained after calibration is discussed in Section 3. Based on published data [9], the specific storage coefficients of the AHU and aquitard were chosen to be constant at 0.003 m⁻¹ and 0.01 m⁻¹, respectively. The specific storage coefficient for the SHU was also part of the calibration of the transient flow model presented in Section 3. Based on the literature data [9], the kinematic porosities for the SHU, aquitard, and AHU were set to 0.3, 0.04, and 0.3, respectively.

2.2.3. Initial and Boundary Conditions of the Numerical Groundwater Model

Flow Model

For the simulation of the transient groundwater flow from 2012 to 2020, we used the map of hydraulic heads computed at steady-state conditions, representing the field situation as observed in 2011. The boundary conditions chosen for flow model 2011 were as follows (Figure 5):

• Constant hydraulic heads (h = 0 m) prescribed at the western boundary (seaside).
• No flow conditions (Q = 0) were applied at the northern boundary, as the selected boundary was approximately perpendicular to the isolines of the hydraulic heads observed in earlier studies.
• No flow conditions (Q = 0) in the south of the western part of the model boundary, as the selected boundary was approximately perpendicular to the isolines of the hydraulic heads observed in earlier studies.
• Constant hydraulic heads were applied to the eastern part of the model boundary in the south, as the chosen boundary limit was not perpendicular to the isolines of the hydraulic heads observed in previous studies.
• Groundwater recharge from rainfall (R [L T⁻¹]) was measured at the meteorological station for 2011 to 2020. The effective aquifer recharge (AR [L T⁻¹]) was calculated considering the potential infiltration coefficient (PIC) [44] for each geological formation in the model domain as follows:
  $${AR}_i=R_i\times {PIC}_i$$

  (3)
  In this study, five specific zones were considered to obtain a local recharge of the groundwater due to rainfall (see Figure 5). For the steady-state flow model, the annual average rainfall and the local PIC values were used to quantify the effective aquifer recharge of each of the five specific zones. For the transient flow model, from 2012 to 2020, the monthly rainfall averages were used. The monthly rainfall averages were then multiplied by the PIC value to obtain the zonal groundwater recharge expressed in m d⁻¹ (see the Supplementary Materials, Figure S1).
• Lateral groundwater recharge (LR) on the eastern boundary: the lateral groundwater recharge was evaluated from the water flux based on the amount of rainfall in the neighboring watersheds, their areas, and their potential infiltration coefficients, as follows:
  $$L{R}_i=C\times {\displaystyle \frac{A_i}{l_i\times d_i}}\times {R_i\times PIC}_i$$

  (4)

  where A_i [L²] is the area of the watershed, l_i [L] and d_i [L] are the length and the average depth of the contact zone between the watershed and model domain, and C is the percentage of lateral groundwater recharge from the total recharge of the watersheds. In our study, coefficient C was determined to be 0.25, as the main amount of rainfall infiltrates vertically in deeper geological formations. The chosen neighboring watersheds are shown in Figure 3.
  In the case of steady-state flow modeling, the lateral groundwater recharge was quantified based on the annual rainfall average of each of the 12 watersheds. In the case where the contact zone between the watershed and model domain presented very different depths, the contact zone was divided into subzones. This resulted in a total of 18 subzones (see the Supplementary Materials, Figure S2). For the transient flow model (2012–2020), the same modeling strategy was used. Based on monthly rainfall averages, the numerical model was built with a time series of incoming flow rates (m d⁻¹) for each of the 18 subzones (see the Supplementary Materials, Figure S3).
• The drainage network was modelled as the inner boundaries on the nodes located nearest to the main ditches, where the hydraulic head was prescribed with a specific water flux constraint (Q < 0). In this flow model, a fixed-head boundary condition was set on the first slice of the model to each node of the main drainage network. In our case, the drainage network represents open ditches. Based on field observations, the hydraulic head (h [L]) values at the defined nodes were considered equal to the elevation of each node (z [L]) at the soil surface minus the estimated depth of the water level of 0.8 m,
  $$h_i=z_i-0.8\mathrm{m} , if Q_i\le Q_{max}=0 {\mathrm{m}}^3{\mathrm{d}}^{-1}$$

  (5)
  The constant-hydraulic-head boundary condition in each of these nodes located along the ditches was therefore constrained by a maximum water flux of 0 m³ d⁻¹, which means that the head boundary condition is only active when water flows out. Hydraulic head conditions above the water level would start infiltrating in cases where the fixed heads are higher than the hydraulic head at the surrounding nodes. Therefore, the constrained-flux condition inhibits inflow.
• A total of 115 wells for industrial use and agricultural use were implemented in the numerical flow model. The two wells of the agricultural and milk processing cooperative were considered, each with a corresponding constant daily extraction rate of 750 and 1400 m³, respectively. Assuming a daily withdrawal of approximately 80 to 100 L of water per dairy cattle, 113 local withdrawal wells for agricultural use were further simulated, with constant flow rates ranging from 6.63 to 38.69 m³ d⁻¹, corresponding to a total flow rate of approximately 2400 m³ d⁻¹.

Transport Model

To simulate the fate of nitrates in the NVZ from 2012 to 2020, the initial conditions for the transport model adopted the observed concentration data for October 2011, which were processed through the interpolation and extrapolation of measured values.

The boundary conditions chosen for the transient transport model were as follows:

• Advective outflow conditions at the western boundary (seaside);
• Impervious conditions for total (convective and dispersive) fluxes applied at the southern and northern boundaries of the model;
• Constant zero concentrations of nitrate were imposed at the eastern boundary nodes of the model domain;
• The total influx of contaminants at the groundwater table of the NVZ (Figure 6) were predefined using the divergence form of the transport equation. Four nitrate input scenarios were considered: (1) scenario N0 corresponds to the initial nitrate concentrations observed in 2011 with no further input (2012–2020); (2) scenario N1 used nitrate input fluxes that corresponded to minimum values of 42 kg N ha⁻¹ yr⁻¹ as observed by [35]; (3) scenario N2 used a nitrate input flux that corresponded to maximum values of 110 kg N ha⁻¹ yr⁻¹, as observed by [35], and (4) scenario N3 was based on the averaged nitrate input fluxes of scenarios N1 and N2 that corresponded to 84 kg N ha⁻¹ yr⁻¹. The masses of N were first converted into NO₃ masses and then into daily NO₃ mass fluxes per square meter. In this study, two periods of fertilization corresponding to the maize and ryegrass culture periods were considered. The period from June to September was expected to account for 32% of the N mass spread (corresponding to a low input of NO₃ at the water table, see Figure 6), and the period from October to January was expected to account for 68% [35].

Macrodispersion in FEFLOW is handled by default by a linear Fickian relationship, distinguishing between longitudinal dispersivity (in the flow direction) and transverse dispersivity (perpendicular to the flow direction). In our transport simulations, we assumed a constant ratio of longitudinal dispersivity to transverse dispersivity of 10, which is a common modeling approach for field-scale applications [9,12]. The longitudinal dispersivity at the field scale, or the longitudinal macrodispersivity, can range from tens to hundreds of meters due to averaging processes in a heterogeneous flow field [9]. Since the triangular grid chosen for the spatial discretization of the model domain contained maximum length sizes of at least 100 to 140 m, we chose a longitudinal macrodispersivity of half the length size of 140 m as the default value. Longitudinal and transverse macrodispersivities of 70 m and 7 m were therefore selected as representative parameters to take into account the effect of field heterogeneities that are not implicitly considered in the spatial distribution of hydraulic permeability.

2.2.4. Numerical Parameters

The dimensionless error criterion in FEFLOW was used for both the iterations in the flow and transport simulations and the automatic time-stepping process. The chosen error tolerance was 10⁻⁴, using the Euclidian L2 integral (RMS) norm. To solve both the symmetric matrix of the fluid-flow problem and the unsymmetric matrix of the nitrate transport problem, the parallel sparse direct solver, PARDISO, implemented in FEFLOW was used.

In the transport model, the fully upwind method was used to improve the numerical stability and reduce the risk of nonphysical negative concentrations. However, this added artificial dispersion to the model as the steep concentration gradients appearing especially at the front of contaminant plumes were smoothed out.

2.3. Climate Change Scenario Analysis

The groundwater flow model was run for both historical input data and data from regional climate models (RCMs) developed as part of the EURO–CORDEX project [45] in collaboration with the University of Parma.

The precipitation and temperature data recorded at Oristano Santa Lucia weather station over the period 1976–2005 were used to correct the systematic errors (bias) in the outputs of the regional climate models. The data showed no gaps in the time series.

RCMs are advantageous for understanding the local climate of regions with complex topography. Historical data were used for downscaling and projection into the future for two RCP scenarios: 4.5 and 8.5. The grid resolution of the RCMs used was 12.5 km (EUR–11 grid), and they were a combination of several general circulation models (GCMs) and RCMs.

Six climate scenarios (out of the seventeen generated) were selected as the most important ones. Each of the six climate scenarios was based on the representative concentration pathways (RCPs) of 4.5 and 8.5, resulting thus in twelve case studies. The regional climate models used for the climate change scenarios are listed in

Table 2. List of RCMs that were used for the climate change scenario analysis.

Model Scenario	Regional Climate Model
1	‘CNRM_CERFACS_CNRM_CM5_CCLM4_8_17’
4	‘DMI_HIRHAM5_NorESM1-M’
8	‘ICHEC_EC_EARTH_HIRHAM5’
9	‘IPSL-INERIS_WRF381P_IPSL-CM5A-MR’
12	‘KNMI_CNRM-CM5’
17	‘MPI_M_MPI_ESM_LR_RCA4’

in order of priority.

Based on the predicted daily rainfall data, the monthly averaged values that were used for the quantification of the effective aquifer recharge using Equation (3) were calculated. Detailed information on the temporal variation of each of the six climate scenarios and for RCPs 4.5 and 8.5 can be found in the Supplementary Materials (Figures S4 and S5).

3. Results and Discussion

3.1. Steady-State Flow Model

A steady-state flow simulation for 2011 was conducted. More than 100 observation wells were used (most of them from the piezometer network of the Regional Environmental Protection Agency, ARPAS) to assess the spatial variation in the hydraulic conductivity of the sandy unit. The observation wells of the ARPAS network are numbered with the suffix “ar”, while the other observation wells belonging to a non-public network, which were used for the calibration of the flow model, are numbered with the suffix “ac”. A trial-and-error method was used to calibrate the numerical model, adjusting the hydraulic conductivities on preselected zones (Figure 7). The good quality of the calibration is shown with the scatter plot of the computed hydraulic heads versus the observed water heads, indicating a low mean error $\overline{\left(E\right)}$, root mean square error (RMS), and standard deviation (σ) of 1.99, 2.98, and 2.99 m, respectively.

Table 3. Water balance components of the steady-state flow model (2011).

		Flow Rates of the Steady-State Flow Model (2011)
		(×103 m3 d−1)
Fixed-head boundary	inflow (+)	5.71
	outflow (−)	128.27
Wells	inflow (+)	-
	outflow (−)	4.55
Lateral groundwater recharge (+)		15.40
Groundwater recharge (+)		111.71

summarizes the individual components of the water balance that were calculated using the calibrated flow model. The groundwater balance computed for 2011 was strongly influenced by the open ditches in the model area. Of the total of approximately 130,000 m³ d⁻¹ of outflowing groundwater, approximately 57% fell on the near-surface drainage system.

This withdrawal rate corresponded to approximately 27 hm³ yr⁻¹, which was approximately 35% of the average withdrawal rates recorded at the Sassu and Luri pumping stations. The local rates of extraction by pumps (agriculture and industry) were only approximately 3.5%. Groundwater was largely fed by natural groundwater recharge via the unsaturated soil zone by up to 85%. A total of 11% of the total inflow came from lateral inflow over the eastern model boundary from the mountain watersheds, and only 4 to 5% of the total inflow rate was from the southern fixed-head boundary.

3.2. Transient Flow Model

In the transient groundwater flow simulation for the period of from 2012 to 2020, the hydraulic heads from the calibrated steady-state groundwater flow model were used as the initial groundwater levels, and the storage coefficients of the different aquifer units were calibrated (Figure 8). The transient boundary conditions were related to the lateral recharge on the eastern boundary and groundwater recharge due to rainfall.

The quality of the calibration is shown in Figure 8 in the form of the scatter plot obtained for more than 100 observation wells located for most of the SHU and AHU units. The hydraulic heads computed at these wells were compared to the groundwater levels measured by ARPAS (Regional Environmental Protection Agency) and by other studies [36,46]. At the end of 2015, the calibration was qualified by a rather low mean error $\overline{\left(E\right)}$, root mean square error (RMS), and standard deviation (σ) of 1.33, 1.71, and 1.72 m, respectively. At the end of the transient flow simulation, the quantified errors were slightly increased by approximately 0.3 to 0.4 m.

Figure 9 shows the computed water heads at four selected observation wells of the ARPAS network placed in the NVZ. For example, for P40ar, the observed time series was rather well represented by the numerical model; P20ar or P10ar were systematically underestimated by approximately half a meter. However, the computed water heads of P46ar were somewhere in between the observed water levels.

The four selected observation wells showed typical results that were also found at the other measuring points as a whole; the computed hydraulic heads showed smaller amplitude fluctuations or larger amplitude fluctuations compared to the measured hydrographs. An offset was visible between the computed and measured hydraulic heads at some measuring points, such as P10ar or P20ar. Due to the lack of small-scale constraints on the lithology properties (possible clay lenses, layered aquifer with different layer widths) and the rough-scale calibration of specific storage, the computed values did not take into account the effect of local heterogeneities, which may have caused different hydraulic head behaviors at the observation wells. An indication of this could be, among other things, the strong trimestral variations in the groundwater readings during the 14 years at monitoring well P20ar (see the Supplementary Materials, Figure S6), which were characterized by a mean hydraulic head and standard deviation of 2.86 m and 0.77 m, respectively.

The influence of the drainage network on the extraction of the near-surface groundwater is shown in Figure 10. The numerical flow model showed, as in the steady-state case, that approximately 58% of the total withdrawal occurred due to drainage (Figure 10a). The total volume of water entering the model area via the southern fixed-head boundaries, even in the transient flow case, was 18 hm³ over the 9 years, representing only approximately 4% of the total water balance. However, the 420 m³ of water computed to have flowed into the model area via near-surface drainage by the model is a numerical artefact. In addition, the results of the numerical flow model showed that the withdrawal rates were relatively constant over the 9 years, with an average of approximately 80,000 m³ d⁻¹, with a slight increase as in the case of the total discharge rates across the fixed-head boundaries (Figure 10b).

3.3. Fate of Nitrates from 2012 to 2020

Considering the monitoring data of the groundwater quality in the NVZ carried out by ARPAS (Regional Environmental Protection Agency) since 2007, as the observed nitrite concentrations and nitrite mass were very low, denitrification and transformation of nitrate in the groundwater can be neglected. In these transport modeling studies, nitrate was therefore assumed to act as a nonreactive solute. Furthermore, the chosen modeling scenarios only comprised a nitrate input at the water table of the NVZ. As only a few observation wells are located in the eastern part of the study area, no input of nitrates into these areas was applied. Furthermore, the detailed locations of the two geological units (the SHU and AHU) are not yet well identified in the southern part of the model domain; the computed mass of NO₃ was therefore only compared to the observed stock of NO₃ in the northern part of the NVZ.

The numerical simulation of the transport of nitrates in the Arborea plain focused on the time period between 2012 and 2020.

The groundwater pollution by nitrates observed at the end of 2011 was used as the initial conditions for the nitrate distribution in the model (Figure 11).

The computed nitrate map shown in Figure 11 is the one obtained with the scenario of zero input of nitrate (at the water table) during the considered period. It is apparent that only parts of the NVZ turned from a “red” color to a “blue” color. There were still local high concentrations of nitrate that reached 150 mg L⁻¹ or higher. As seen from the scatter plot shown in Figure 12, a zero-nitrate input at the water table contributed to computed concentrations that largely underestimated the nitrate concentrations observed in the ARPAS observation wells in December 2015 and December 2020. Furthermore, the quantified root mean square error (RMS) increased significantly from 73.29 to 85.53 mg L⁻¹ from the end of December 2015 to 2020.

Figure 13 shows the nitrate concentrations computed based on the chosen scenarios at four selected observation wells placed in the NVZ. Modeling scenario N2, corresponding to the maximum nitrate input, came closest to the observations but still underestimated the real stock of nitrate in the sandy unit of the NVZ estimated at the end of December 2020 (

Table 4. Observed and computed stock of nitrate mass in the SHU of the northern part of the NVZ.

2011—NO3 Mass Observed [×103 kg]	End of 2020—NO3 Mass Observed [×103 kg]	End of 2020—NO3 Mass Computed [×103 kg]
11,100	13,915	Zero input (SimN0):	7500
		Min input (SimN1):	9128
		Max input (SimN2):	11,110
		Avg input (SimN3):	10,335

). The differences between the calculated time-dependent concentrations and the measured concentrations at the four selected monitoring wells were most likely due to the spatio-temporal variation in nitrate inputs at the soil surface within the NVZ, the transition time for the dissolved nitrate to reach the water table, and the unknown (and not yet accounted for) nitrate input to the groundwater east of the NVZ. Furthermore, as an initial modeling approach, our simulations of the fate of nitrate in the aquifer were based on data from field experiments on locally observed nitrate fluxes and did not aim to obtain an adjusted map of the spatial variation in the nitrate input at the water table as a final result. This could be part of future field investigations and modeling studies.

3.4. Modeling of Climate Change Scenarios: Groundwater Levels in the near Future, Medium Term, and Long Term and the Fate of Nitrate Pollution in the near Future

The time period of the modeled scenarios was 2021–2098. To implement the selected 12 scenarios (see Table 2), time-varying monthly averaged groundwater recharge values in five areas were used, which were derived from the predicted rainfall of the RGM–GCM models for both RCPs 4.5 and 8.5. The pumping rates of the extraction wells and the drainage network were maintained as per the transient flow model for the simulation period (2012–2020). It is worth noting that for reasons of simplification, the spatially variable distribution of the groundwater inflow at the eastern model boundary was maintained, with the constant inflow rates calculated as the time averages from 2012 to 2020.

Table 5. Mean annual rainfall at Oristano Santa Lucia weather station as a function of time and RCP considered.

	Near-Future (2021–2040) Rainfall	Medium-Term (2041–2060) Rainfall	Long-Term (2081–2098) Rainfall
Mean annual rainfall (mm) for RCP 4.5	574.07	574.70	604.15
Mean annual rainfall (mm) for RCP 8.5	566.57	579.39	545.38

shows the mean annual rainfall predicted at Oristano Santa Lucia weather station as a function of the period under consideration and the RCP considered. Until 2060, the mean annual rainfall predicted for RCP 4.5 remained almost constant at 574 mm and was fairly close to the value predicted for RCP 8.5. In the long-term perspective (2081–2098), however, the annual rainfall increased significantly for RCP 4.5 and decreased for RCP 8.5.

The influence of climate change scenario 1 on the groundwater level of four of the selected monitoring wells is exemplarily shown in Figure 14. High groundwater fluctuations were only be seen in the area of the upward-coning (low-permeability) aquitard (see P40ar). Otherwise, the hydraulic heads computed in the NVZ did not fluctuate greatly over time, as the groundwater table is constrained by the dense network of artificial ditches.

More detailed information on the temporal variation in the hydraulic heads computed at the individual ARPAS measuring points can be found in the Supplementary Materials. Tables S1 and S2 summarize the temporal averages of the minimum, mean, and maximum values of the hydraulic heads calculated for the three time periods and for the 43 monitoring wells based on the six scenarios for RCP 4.5 and RCP 8.5, respectively. The three time periods considered were the near future (2021–2040), medium-term perspective (2041–2060), and long-term perspective (2081–2098).

Figure 15 summarizes an interesting result regarding the long-term perspective (2081–2098). The deviation in the mean hydraulic head was determined from the six scenarios of the hydraulic head computed in December 2020 for RCPs 4.5 and 8.5. Based on the 43 monitoring wells placed in the NVZ, the hydraulic heads for RCP 8.5 decreased significantly (by several tens of cm) compared to RCP 4.5. This was probably due to the fact that the average annual rainfall increased by 6% for RCP 4.5 and decreased by 4% for RCP 8.5 (see Table 5).

To analyze the future evolution of nitrate pollution in the near future (end of 2040), initial transport modeling was performed using the scenario 1 for RCP 4.5 and applying zero nitrate at the NVZ water table (Figure 16). The hotspots with nitrate concentrations higher than 200 mg L⁻¹ observed at the end of 2020 were slightly attenuated 20 years later, but concentrations of 50 mg L⁻¹ and more were still present in the Arborea coastal aquifer.

Figure 17 shows the computed nitrate concentrations at 10 observation wells in the NVZ, where the concentrations at the end of the near-future scenario (December 2040) were between 50 and 100 mg L⁻¹.

At observation well P04ar, where very high nitrate concentrations of up to 400 mg L⁻¹ were measured between 2012 and 2017, the zero-nitrate scenario and the advective–dispersive effect of the groundwater flow enabled a significant reduction in the concentration from 200 mg L⁻¹ in December 2020 to 100 mg L⁻¹ at the end of the simulation period. Although a reduction in the nitrate concentration was achieved at most of the monitoring wells, a slight increase in the nitrate concentration up to 50 mg L⁻¹ was recorded at some of the monitoring wells, such as P03ar. This was most likely due to the shift in the higher nitrate concentrations from upstream.

4. Conclusions

The developed partially unstructured transient flow model (2012–2020) took into account the seasonal recharge of the aquifer system by rainfall and lateral inflow at the model boundary, the drainage network, and the water volumes pumped for industrial use and by farms for technical use and livestock. In these flow simulations, local potential infiltration coefficients, which depend on the soil surface texture multiplied by the monthly precipitation, were used to assess both the effective vertical and lateral transient groundwater recharge. Compared to the full hydrologic water balance approach used in other field studies, this modeling approach is much simpler in practical applications. For example, to evaluate the vertical groundwater recharge, Matzeu [39] derived the transient variation in the vertical inflow rate at the groundwater table from a water balance analysis of the unsaturated zone. Furthermore, our modeling approach can be easily implemented in existing groundwater flow models and effortlessly allows for the direct integration of monthly rainfall data from climate change scenarios to predict the future evolution of groundwater reserves.

The calibrated steady-state flow model representative of 2011 indicated that groundwater recharge due to rainfall was the main component of the water input. The lateral recharge at the eastern boundary corresponded to approximately 13% of the annual water budget, whereas the industrial wells together with the wells used by farmers corresponded to only 3.5%. The results of the transient flow model highlighted the influence of the drainage network on the overall groundwater management; nearly steady-state flow rates were predicted by the numerical model, with the total water volume drained by ditches accounting for approximately 58% of the annual outflow volume. It is worth noting that this rather high withdrawal rate is of the order of magnitude of the average withdrawal rates recorded at the Sassu and Luri pumping stations.

The numerical transport simulations conducted from 2012 to 2020 using hypothetical field-based nitrate input scenarios globally underestimated the high concentrations observed in the NVZ. The reason for this is most likely the working hypothesis that the current simulation scenarios neglected the application of NO₃ in the eastern part of the model area. However, as observed in the field, the computed nitrate concentrations in December 2020 still varied strongly in space, from several mg L⁻¹ to several hundreds of mg L⁻¹. The origin of these local hotspots, whose computed concentrations are certainly somewhat depleted, is not yet known. The final computed concentration field had a pattern similar to the one observed in 2011. Furthermore, the total mass of nitrates remaining in the aquifer observed in December 2020 was approximately 20% higher than the stock of nitrates estimated in 2011. The reason for this general trend is still uncertain, but one possibility might be the storage of nitrate in clay lenses and small layered inclusions that were not explicitly taken into account in the numerical model domain.

According to the numerical groundwater model, it seems that denitrification and transformation of nitrate can be neglected, as the observed nitrite concentrations and nitrite mass were very low. However, previous studies using isotopic methods have noted that locally significant denitrification processes exist.

It is worth noting that groundwater level fluctuations based on the climate change scenario and zero-nitrate input scenario reduced the nitrate pollution; hotspots were attenuated, but concentrations in several monitoring wells still exceeded the limit of 50 mg L⁻¹.

Ongoing research is aimed at quantifying the spatiotemporal distribution of nitrate in the SHU aquifer under transient groundwater flow conditions to compare different water management methods, climate change scenarios, and contamination scenarios. As recently discussed with stakeholders, these numerical studies could be combined with field studies aiming to test some “upstream” and “downstream” mitigation options to reduce the nitrate groundwater contamination in the Arborea plain. “Upstream” options include modifications to livestock, fertilization, and crop management in order to reduce the input of nitrate in the soil and aquifers, while “downstream” options refer to methods to restore the polluted aquifer, such the Forested Infiltration Area (FIA) technique, by exploiting, for example, the potential positive effect of eucalyptus hedges on groundwater denitrification (the MENAWARA, NATMed, and SARNITRO projects) [47,48].

This study also highlights that in such an intensively used coastal aquifer system, the collection of detailed local data regarding groundwater level fluctuations and possible small-scale aquifer heterogeneities and existing denitrification conditions is of paramount importance to allow for an improved understanding of the overall system and, thus, the establishment of an appropriate water management plan.