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

Gerhard Schäfer, Manon Lincker, Antonio Sessini, Alberto Carletti

► To cite this version:

Gerhard Schäfer, Manon Lincker, Antonio Sessini, Alberto Carletti. Numerical Groundwater Model to Assess the Fate of Nitrates in the Coastal Aquifer of Arborea (Sardinia, Italy). Water, 2024, 16 (19), pp.2729. 10.3390/w16192729. [https://hal-05407065]

HAL Id: hal-05407065

https://hal.science/hal-05407065v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

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

Gerhard Schäfer<fnref n="1" />, Manon Lincker<fnref n="1" />, Antonio Sessini<fnref n="2" /> and Alberto Carletti<fnref n="3" />,*

**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⁻<fnref n="1" /> to several hundreds of mg L⁻<fnref n="1" />. 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⁻<fnref n="1" />.

**Keywords:** groundwater recharge/water budget; drainage network; nitrate pollution; climate change; Italy

## 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

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Schäfer, G.; Lincker, M.; Sessini, A.; Carletti, A. Numerical Groundwater Model to Assess the Fate of Nitrates in the Coastal Aquifer of Arborea (Sardinia, Italy). *Water* **2024**, *16*, 2729. https://doi.org/10.3390/w16192729

Academic Editor: Andrzej Witkowski

Received: 29 July 2024

Revised: 21 September 2024

Accepted: 23 September 2024

Published: 25 September 2024

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Institut Terre et Environnement de Strasbourg (ITES), UMR 7063 CNRS-Université de Strasbourg, 67000 Strasbourg, France; schafer@unistra.fr (G.S.); manon.lincker@live.com (M.L.)</footnote>

<footnote n="2">² Department of Chemical and Geological Sciences, University of Cagliari, Cittadella Universitaria di Monserrato, Blocco A—S.P. Monserrato-Sestu, km 0.700, 09042 Monserrato, Italy; antoniosessini@tiscali.it</footnote>

<footnote n="3">³ Department of Civil, Environmental Engineering and Architecture, University of Cagliari, Via Marengo 2, 09123 Cagliari, Italy</footnote>

<footnote>* Correspondence: alberto.carletti@unica.it</footnote>
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–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–18]. Furthermore, the soil lithology, soil texture, and, particularly, the depth of the water table play critical roles in groundwater recharge [19–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–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–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 Pani-
coni [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 evalua-
tion 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 Marceddi 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].

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Location of the study area. The municipal boundaries of Arborea are shown in yellow, and the boundaries of the nitrate vulnerable zone (NVZ) perimeter are shown in red.</figcaption></figure>

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^{-5}$ and $10^{-6}$ m s$^{-1}$. The AHU is a multilayer aquifer with good permeability ($K = 10^{-4}-10^{-5}$ m s$^{-1}$) 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.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Hydrogeology of the study area in 3D (modified from Ghiglieri et al. [36]). Two main aquifers are identified; 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 an aquitard of silt and clay that outcrops in the Sassu Lagoon. The aquitard has been shown to become thinner towards the south of the area until both aquifers come into contact.</figcaption></figure>

### 2.1.3. Water Use

The Arborea farming district receives high-pressure irrigation water from the Eleonora d'Arborea dam (some 35 $hm^3 yr^{-1}$), located approximately 35 km away. The Oristanese Land Reclamation Consortium distributes an annual amount of surface water equal to 28.4 $hm^3 yr^{-1}$ at the field scale by a permanent sprinkler system for a total irrigated area of 4900 ha (approximately 5800 $m^3 ha^{-1}$). Groundwater, mainly from the SHU aquifer, is used by farmers only for cattle watering (approximately 0.97 $hm^3 yr^{-1}$) and for the washing of livestock facilities (approximately 0.04 $hm^3 yr^{-1}$) (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^3 yr^{-1}$ and 0.27 $hm^3 yr^{-1}$, 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, 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^{-1}$. Regarding the other N compounds, only 4 sampling points (9%) showed average annual ammonium concentrations above the threshold value of 0.5 $mg L^{-1}$, while 14 of them (32%) had average annual nitrite concentrations up to 6 times higher than the threshold limit (0.5 $mg L^{-1}$).
<table><caption>**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.</caption><thead><tr><th rowspan="3">ID_WELL<br>(ARPAS Network)</th><th rowspan="3">X Coordinate<br>(WGS 84)</th><th rowspan="3">Y Coordinate<br>(WGS 84)</th><th colspan="2">NO<sub>3</sub></th><th colspan="2">NH<sub>4</sub></th><th colspan="2">NO<sub>2</sub></th></tr><tr><th>Average</th><th>SD</th><th>Average</th><th>SD</th><th>Average</th><th>SD</th></tr><tr><th>mg L<sup>-1</sup></th><th>mg L<sup>-1</sup></th><th>mg L<sup>-1</sup></th><th>mg L<sup>-1</sup></th><th>mg L<sup>-1</sup></th><th>mg L<sup>-1</sup></th><th>mg L<sup>-1</sup></th></tr></thead><tbody><tr><td>17PZ008ar</td><td>461,524</td><td>4,407,087</td><td>11.55</td><td>11.78</td><td><strong>0.75</strong></td><td>0.64</td><td>0.02</td><td>0.00</td></tr><tr><td>P01ar</td><td>461,793</td><td>4,404,058</td><td>12.88</td><td>7.80</td><td>0.08</td><td>0.06</td><td>0.43</td><td>0.53</td></tr><tr><td>P02ar</td><td>462,647</td><td>4,404,165</td><td><strong>171.51</strong></td><td>26.79</td><td>0.12</td><td>0.12</td><td><strong>1.38</strong></td><td>0.85</td></tr><tr><td>P03ar</td><td>463,755</td><td>4,404,111</td><td>23.76</td><td>26.02</td><td>0.26</td><td>0.24</td><td>0.30</td><td>0.83</td></tr><tr><td>P04ar</td><td>461,798</td><td>4,402,316</td><td><strong>329.81</strong></td><td>59.81</td><td>0.17</td><td>0.25</td><td><strong>3.06</strong></td><td>1.44</td></tr><tr><td>P05ar</td><td>461,060</td><td>4,401,284</td><td>32.93</td><td>38.13</td><td>0.12</td><td>0.13</td><td><strong>0.63</strong></td><td>0.48</td></tr><tr><td>P06ar</td><td>462,012</td><td>4,401,286</td><td><strong>98.45</strong></td><td>54.44</td><td>0.34</td><td>0.43</td><td><strong>0.72</strong></td><td>0.45</td></tr><tr><td>P07ar</td><td>462,652</td><td>4,401,289</td><td>30.05</td><td>15.32</td><td>0.08</td><td>0.05</td><td>0.32</td><td>0.21</td></tr><tr><td>P08ar</td><td>463,609</td><td>4,401,302</td><td><strong>97.76</strong></td><td>37.43</td><td>0.05</td><td>0.03</td><td>0.18</td><td>0.10</td></tr><tr><td>P09ar</td><td>464,530</td><td>4,400,935</td><td><strong>83.07</strong></td><td>69.89</td><td><strong>1.02</strong></td><td>1.33</td><td>0.24</td><td>0.17</td></tr><tr><td>P10ar</td><td>461,269</td><td>4,399,687</td><td><strong>125.11</strong></td><td>23.07</td><td>0.16</td><td>0.23</td><td><strong>1.94</strong></td><td>0.78</td></tr><tr><td>P11ar</td><td>463,166</td><td>4,398,687</td><td><strong>115.62</strong></td><td>37.07</td><td>0.13</td><td>0.11</td><td><strong>0.97</strong></td><td>0.74</td></tr><tr><td>P12ar</td><td>462,203</td><td>4,397,650</td><td><strong>189.43</strong></td><td>74.92</td><td>0.06</td><td>0.10</td><td>0.09</td><td>0.13</td></tr><tr><td>P13ar</td><td>462,124</td><td>4,401,291</td><td><strong>76.25</strong></td><td>29.89</td><td>0.12</td><td>0.20</td><td><strong>0.90</strong></td><td>0.46</td></tr><tr><td>P14ar</td><td>462,975</td><td>4,400,705</td><td><strong>181.03</strong></td><td>44.61</td><td>0.04</td><td>0.05</td><td><strong>0.56</strong></td><td>0.42</td></tr><tr><td>P15ar</td><td>464,029</td><td>4,401,664</td><td><strong>89.02</strong></td><td>30.46</td><td>0.06</td><td>0.06</td><td>0.31</td><td>0.37</td></tr><tr><td>P16ar</td><td>464,276</td><td>4,401,706</td><td>43.13</td><td>23.78</td><td>0.06</td><td>0.17</td><td>0.42</td><td>0.27</td></tr><tr><td>P17ar</td><td>462,650</td><td>4,398,853</td><td>44.09</td><td>15.71</td><td>0.06</td><td>0.07</td><td><strong>0.78</strong></td><td>0.60</td></tr><tr><td>P20ar</td><td>463,601</td><td>4,406,282</td><td><strong>286.70</strong></td><td>70.92</td><td>0.07</td><td>0.06</td><td><strong>0.60</strong></td><td>0.62</td></tr><tr><td>P21ar</td><td>459,789</td><td>4,398,871</td><td>22.01</td><td>31.33</td><td>0.32</td><td>0.20</td><td><strong>0.72</strong></td><td>0.75</td></tr><tr><td>P22ar</td><td>461,451</td><td>4,399,488</td><td><strong>148.99</strong></td><td>37.82</td><td>0.14</td><td>0.16</td><td><strong>2.07</strong></td><td>0.84</td></tr><tr><td>P23ar</td><td>458,950</td><td>4,397,389</td><td><strong>198.40</strong></td><td>84.22</td><td>0.12</td><td>0.15</td><td><strong>2.13</strong></td><td>1.30</td></tr><tr><td>P24ar</td><td>461,293</td><td>4,396,566</td><td><strong>57.34</strong></td><td>32.41</td><td><strong>0.77</strong></td><td>1.41</td><td>0.48</td><td>0.30</td></tr><tr><td>P26ar</td><td>464,449</td><td>4,406,075</td><td><strong>226.38</strong></td><td>26.00</td><td>0.03</td><td>0.03</td><td>0.04</td><td>0.02</td></tr><tr><td>P27ar</td><td>464,024</td><td>4,406,258</td><td><strong>64.33</strong></td><td>12.57</td><td>0.05</td><td>0.06</td><td>0.42</td><td>0.17</td></tr><tr><td>P28ar</td><td>463,399</td><td>4,401,307</td><td><strong>80.38</strong></td><td>23.61</td><td>0.06</td><td>0.07</td><td>0.32</td><td>0.08</td></tr><tr><td>P29ar</td><td>462,651</td><td>4,400,515</td><td>2.96</td><td>5.22</td><td>0.21</td><td>0.07</td><td>0.12</td><td>0.34</td></tr><tr><td>P30ar</td><td>460,104</td><td>4,398,210</td><td>3.06</td><td>7.81</td><td><strong>1.18</strong></td><td>0.19</td><td>0.15</td><td>0.25</td></tr><tr><td>P31ar</td><td>462,498</td><td>4,406,741</td><td>2.50</td><td>2.75</td><td>0.21</td><td>0.08</td><td>0.04</td><td>0.03</td></tr><tr><td>P32ar</td><td>462,838</td><td>4,406,745</td><td>1.50</td><td>1.09</td><td>0.15</td><td>0.11</td><td>0.03</td><td>0.02</td></tr><tr><td>P34ar</td><td>461,572</td><td>4,400,003</td><td><strong>53.38</strong></td><td>32.38</td><td>0.07</td><td>0.07</td><td>0.17</td><td>0.17</td></tr><tr><td>P35ar</td><td>462,881</td><td>4,402,472</td><td><strong>68.51</strong></td><td>16.49</td><td>0.08</td><td>0.11</td><td><strong>0.98</strong></td><td>0.71</td></tr><tr><td>P36ar</td><td>460,342</td><td>4,400,943</td><td>40.04</td><td>18.60</td><td>0.08</td><td>0.06</td><td>0.43</td><td>0.27</td></tr><tr><td>P37ar</td><td>460,418</td><td>4,400,804</td><td>5.73</td><td>19.21</td><td>0.19</td><td>0.07</td><td>0.09</td><td>0.09</td></tr><tr><td>P38ar</td><td>460,691</td><td>4,396,128</td><td><strong>141.24</strong></td><td>46.03</td><td>0.07</td><td>0.13</td><td>0.12</td><td>0.06</td></tr><tr><td>P39ar</td><td>465,035</td><td>4,404,262</td><td><strong>56.37</strong></td><td>38.81</td><td>0.04</td><td>0.02</td><td>0.05</td><td>0.03</td></tr><tr><td>P40ar</td><td>464,805</td><td>4,402,660</td><td><strong>53.73</strong></td><td>16.24</td><td>0.03</td><td>0.01</td><td>0.10</td><td>0.05</td></tr><tr><td>P41ar</td><td>464,823</td><td>4,399,386</td><td>32.35</td><td>18.71</td><td>0.10</td><td>0.14</td><td>0.13</td><td>0.27</td></tr><tr><td>P42ar</td><td>461,452</td><td>4,396,954</td><td><strong>56.97</strong></td><td>20.58</td><td>0.03</td><td>0.01</td><td>0.10</td><td>0.08</td></tr><tr><td>P43ar</td><td>459,512</td><td>4,396,949</td><td><strong>84.23</strong></td><td>24.77</td><td>0.38</td><td>0.66</td><td>0.36</td><td>0.24</td></tr><tr><td>P45ar</td><td>461,841</td><td>4,396,172</td><td>3.25</td><td>2.58</td><td>0.29</td><td>0.61</td><td>0.10</td><td>0.14</td></tr><tr><td>P46ar</td><td>463,935</td><td>4,397,919</td><td>5.54</td><td>1.91</td><td>0.03</td><td>0.01</td><td>0.06</td><td>0.04</td></tr><tr><td>P47ar</td><td>461,109</td><td>4,397,773</td><td>14.02</td><td>23.53</td><td>0.15</td><td>0.09</td><td>0.25</td><td>0.34</td></tr></tbody></table>

## 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 \frac{\partial h}{\partial t} = \frac{\partial}{\partial x_i} \left( K_{ij} \frac{\partial h}{\partial x_j} \right) - Q \quad (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]

$$\frac{\partial C}{\partial t} = -\frac{\partial}{\partial x_i}(v_i C) + \frac{\partial}{\partial x_i}\left(D_{ij}\frac{\partial C}{\partial x_j}\right) + D_m \frac{\partial^2 C}{\partial x_i^2} - QC_{in} \\ \text{with } v_i = -\frac{1}{n_c}K_{ij}\frac{\partial h}{\partial x_j}, \quad (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.

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Model domain, its hydrography network, and associated watersheds at the eastern boundary contributing to the lateral recharge of the upper part of the main aquifer. The contour of the model domain is shown in yellow; numbers 1 to 12 correspond to the identifying number of each watershed.</figcaption></figure>

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.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Setup of the numerical groundwater model based on the hydrogeological model developed by Ghiglieri et al. [36]. The 4182 highlighted elements shown at the ground surface represent the outcrop of the aquitard.</figcaption></figure>

### 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):

1. Constant hydraulic heads ($h = 0$ m) prescribed at the western boundary (seaside).

2. 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.

3. 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.

4. 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.

5. 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 \quad (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).

6. 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:

$$LR_i = C \times \frac{A_i}{l_i \times d_i} \times R_i \times PIC_i \quad (4)$$

where Aᵢ [L²] is the area of the watershed, lᵢ [L] and dᵢ [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).

7. 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.8m, \quad if \ Q_i \le Q_{max} = 0 \ m^3 d^{-1} \quad (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.

8. 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⁻¹.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Boundary conditions used for steady-state flow model 2011: (a) drainage network (nodes are highlighted in yellow), (b) prescribed hydraulic heads ($h = 0$) along the sea side (nodes are highlighted in yellow), (c) zonal distributed ground water recharge, and (d) lateral recharge along the eastern boundary (nodes are highlighted in yellow).</figcaption></figure>

## 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:

* Advection 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].
<figure><img src="image_10.png" /><figcaption>**Figure 6.** NO₃ input applied at the water table of the sandy hydrogeological unit for modeling scenarios N1, N2, and N3. The green and pink lines shown in the top left corner of Figure 6 indicate the two periods of three months in which a low and a high nitrate input into the groundwater took place for the first year of transport modeling.</figcaption></figure>

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^{-4}$, 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 in order of priority.

<table><caption>**Table 2.** List of RCMs that were used for the climate change scenario analysis.</caption><thead><tr><td>Model Scenario</td><td>Regional Climate Model</td></tr></thead><tbody><tr><td>1</td><td>'CNRM_CERFACS_CNRM_CM5_CCLM4_8_17'</td></tr><tr><td>4</td><td>'DMI_HIRHAM5_NorESM1-M'</td></tr><tr><td>8</td><td>'ICHEC_EC_EARTH_HIRHAM5'</td></tr><tr><td>9</td><td>'IPSL-INERIS_WRF381P_IPSL-CM5A-MR'</td></tr><tr><td>12</td><td>'KNMI_CNRM-CM5'</td></tr><tr><td>17</td><td>'MPI_M_MPI_ESM_LR_RCA4'</td></tr></tbody></table>

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 (E), root mean square error (RMS), and standard deviation (σ) of 1.99, 2.98, and 2.99 m, respectively.

Table 3 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.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Steady-state flow simulation: (a) map of calibrated hydraulic conductivities of the upper aquifer (note: the purple patch corresponds to the hydraulic conductivity of the outcropping aquitard), (b) computed hydraulic heads, and (c) scatter plot of computed hydraulic heads versus observed water heads. The blue line shown in (c) represents the 1:1 line.</figcaption></figure>

<table><caption>**Table 3.** Water balance components of the steady-state flow model (2011).</caption><thead><tr><th colspan="3">Flow Rates of the Steady-State Flow Model (2011)</th></tr><tr><th></th><th></th><th>(×10<sup>3</sup> m<sup>3</sup> d<sup>-1</sup>)</th></tr></thead><tbody><tr><td rowspan="2">Fixed-head boundary</td><td>inflow (+)</td><td>5.71</td></tr><tr><td>outflow (-)</td><td>128.27</td></tr><tr><td rowspan="2">Wells</td><td>inflow (+)</td><td>-</td></tr><tr><td>outflow (-)</td><td>4.55</td></tr><tr><td>Lateral groundwater recharge (+)</td><td></td><td>15.40</td></tr><tr><td>Groundwater recharge (+)</td><td></td><td>111.71</td></tr></tbody></table>

This withdrawal rate corresponded to approximately 27 $hm^3$ yr$^{-1}$, 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.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Transient flow simulation: (a) map of calibrated storage coefficients, (b) computed hydraulic heads for December 2020, and scatterplot of computed hydraulic heads and observed water heads at (c) the end of 2015 and (d) the end of 2020. The blue lines shown in (c,d) represent the 1:1 lines.</figcaption></figure>

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 ($\bar{E}$), root mean square error (RMS), and standard deviation ($\sigma$) 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.

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Hydraulic heads computed in the sandy hydrogeological unit of the aquifer from 2012 to 2020 compared to those measured at four observation wells located in the NVZ. Note: t = 0 d and t = 3285 d correspond to 1 January 2012 and 31 December 2020, respectively.</figcaption></figure>

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^3$ over the 9 years, representing only approximately 4% of the total water balance. However, the 420 $m^3$ 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^3 d^{-1}$, with a slight increase as in the case of the total discharge rates across the fixed-head boundaries (Figure 10b).
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Effect of drains on the transient water balance: (a) cumulative water volume computed (in +/out −) across the fixed-head boundaries, (b) outflow rates computed over the fixed-head boundaries for the total domain and the drains. Note: t = 0 d and t = 3285 d correspond to 1 January 2012 and 31 December 2020, respectively.</figcaption></figure>

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 NO3 was therefore only compared to the observed stock of NO3 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).
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Contour map of nitrate distribution in the sandy unit: (a) concentrations observed in October 2011 and (b) concentrations computed at the end of December 2020 based on scenario N0.</figcaption></figure>

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><img src="image_16.png" /><figcaption>**Figure 12.** Scatter plot of computed and observed nitrate concentrations: (a) December 2015 and (b) December 2020. The blue lines shown represent the 1:1 lines.</figcaption></figure>
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). 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.

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Computed nitrate concentrations versus measurements at four selected ARPAS observation wells from 2012 to 2020. Note: t = 0 d and t = 3285 d correspond to 1 January 2012 and 31 December 2020, respectively.</figcaption></figure>

<table><caption>**Table 4.** Observed and computed stock of nitrate mass in the SHU of the northern part of the NVZ.</caption><thead><tr><th>2011—NO<sub>3</sub> Mass<br>Observed [×10<sup>3</sup> kg]</th><th>End of 2020—NO<sub>3</sub> Mass<br>Observed [×10<sup>3</sup> kg]</th><th>End of 2020—NO<sub>3</sub> Mass Computed<br>[×10<sup>3</sup> kg]</th></tr></thead><tbody><tr><td>11,100</td><td>13,915</td><td>Zero input (SimN0): 7500<br>Min input (SimN1): 9128<br>Max input (SimN2): 11,110<br>Avg input (SimN3): 10,335</td></tr></tbody></table>
### 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 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.

<table><caption>**Table 5.** Mean annual rainfall at Oristano Santa Lucia weather station as a function of time and RCP considered.</caption><thead><tr><th></th><th>Near-Future<br>(2021–2040)<br>Rainfall</th><th>Medium-Term<br>(2041–2060)<br>Rainfall</th><th>Long-Term<br>(2081–2098)<br>Rainfall</th></tr></thead><tbody><tr><td>Mean annual rainfall (mm) for RCP 4.5</td><td>574.07</td><td>574.70</td><td>604.15</td></tr><tr><td>Mean annual rainfall (mm) for RCP 8.5</td><td>566.57</td><td>579.39</td><td>545.38</td></tr></tbody></table>

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 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⁻¹.
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Time-dependent hydraulic heads computed in scenario 1 for RCPs 4.5 and 8.5 at four selected monitoring wells of the NVZ from January 2021 to December 2098.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 15.** Histogram of hydraulic head differences calculated for RCPs 4.5 and 8.5 in the long-term perspective.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Figure 16.** Contour map of nitrate concentrations computed for January 2012, December 2020, and December 2040. Nitrate concentrations exceeding 100 mg L⁻¹ are shown in red.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 17.** Computed nitrate concentrations at 10 selected observation points of the NVZ from January 2012 to December 2040.</figcaption></figure>

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 ground-
water 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 hypothet-
ical 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 NO3 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-1 to several hundreds of
mg L-1. 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 as 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.

**Supplementary Materials:** The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/w16192729/s1, Figure S1: Monthly zonal groundwater recharges expressed in m d⁻¹ used for the transient flow model from 2012 to 2020. Note: t = 0 d and t = 3285 d correspond to 1 January 2012 and 31 December 2020, respectively; Figure S2: Steady-state flow model: lateral groundwater recharge (m d⁻¹) used for the chosen 18 subzones; Figure S3: Transient flow model: time series of monthly lateral groundwater recharge (m d⁻¹) used for the chosen 18 subzones. Note: t = 0 d and t = 3285 d correspond to 1 January 2012 and 31 December 2020, respectively; Figure S4: Monthly averaged rainfall data predicted for scenarios 1, 4, 8, 9, 12, and 17 for RCP 4.5 from January 2021 to December 2098; Figure S5: Monthly averaged rainfall data predicted for scenarios 1, 4, 8, 9, 12, and 17 for RCP 8.5 from January 2021 to December 2098; Figure S6: Trimestral groundwater level readings at monitoring well P20ar from 2007 to 2020. The dashed line indicates the observed trend, showing slightly decreasing hydraulic heads; Table S1: Results of selected scenarios 1, 4, 8, 9, 12, and 17 (see Table 2) based on RCP 4.5: time averages of the minimum, mean, and maximum values of computed hydraulic heads; Table S2: Results of selected scenarios 1, 4, 8, 9, 12, and 17 (see Table 2) based on RCP 8.5: time averages of the minimum, mean, and maximum values of computed hydraulic heads.

**Author Contributions:** Conceptualization, A.C., M.L. and G.S.; data curation, A.C., A.S., M.L. and G.S.; modeling, M.L., A.S. and G.S.; formal analysis, M.L. and G.S.; funding acquisition, A.C. and G.S.; methodology, G.S.; project administration, A.C. and G.S.; supervision, A.C. and G.S.; visualization, M.L. and G.S.; writing—original draft, G.S.; writing—review and editing, A.C., M.L., A.S. and G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research study was funded by the General Secretariat for Research and Technology of the Ministry of Development and Investments under the PRIMA program. PRIMA is supported by the Art.185 initiative and co-funded under Horizon 2020, the European Union's Program for Research and Innovation. We also acknowledge funding from the Italian Ministry of University and Research, CUP no. J84D18000180005.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Groundwater levels, rainfall data, and groundwater quality data were collected by ARPAS (Agenzia Regionale per la Protezione dell'Ambiente della Sardegna) and made available for the study, for which the authors express their sincere gratitude. We extend our gratitude to Yves Armando of ITES for analyzing, filtering, and numerically implementing the rainfall data on climate change. The authors would also like to dedicate this study to the memory of their colleague, Professor Giorgio Ghiglieri, who prematurely and suddenly passed away in August 2021 and who was instrumental in the initial scientific steps of the work.

**Conflicts of Interest:** The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## References

1. MED Groundwater Working Group. Mediterranean Groundwater Report. 2007. Available online: https://circabc.europa.eu/sd/a/50c3b2a9-4816-4ab1-9a33-d41c327759e3/Mediterranean%20Groundwater%20Report_final_150207_clear.pdf (accessed on 26 January 2023).

2. UNEP-MAP; UNESCO-IHP. Final Report on Mediterranean Coastal Aquifers and Groundwater Including the Coastal Aquifer Supplement to the TDA-MED and the Sub-Regional Action Plans. Paris: Strategic Partnership for the Mediterranean Sea Large Marine Ecosystem (MedPartnership). 2015. Available online: https://unesdoc.unesco.org/images/0023/002353/235306e.pdf (accessed on 26 January 2023).

3. Rasmussen, P.; Sonnenborg, T.O.; Goncear, G.; Hinsby, K. Assessing impacts of climate change. sea level rise. and drainage canals on saltwater intrusion to coastal aquifer. *Hydrol. Earth Syst. Sci.* **2013**, *17*, 421–443. [CrossRef]

4. Lejeusne, C.; Chevaldonné, P.; Pergent-Martini, C.; Boudouresque, C.F.; Pérez, T. Climate change effects on a miniature ocean: The highly diverse, highly impacted Mediterranean Sea. *Trends Ecol. Evol.* **2010**, *25*, 250–260. [CrossRef] [PubMed]
5. Anderson, M.P.; Woessner, W.W.; Hunt, R.J. *Applied Groundwater Modeling: Simulation of Flow and Advective Transport*, 2nd ed.; Academic Press: London, UK, 2015; pp. 3–25.

6. Ebraheem, A.M.; Riadn, R.; Wycisk, P.; Sefelnasr, A.M. A local-scale groundwater flow model for groundwater resources management in Dakhla Oasis, SW Egypt. *Hydrogeol. J.* **2004**, *12*, 714–722. [CrossRef]

7. Psarropoulou, E.; Karatzas, G. Transient groundwater modelling with spatio-temporally variable fluxes in a complex aquifer system: New approach in defining boundary conditions for a transient flow model. *Civ. Eng. Environ. Syst.* **2012**, *29*, 1–21. [CrossRef]

8. Bear, J. *Hydraulics of Groundwater*; McGraw-Hill Book: New York, NY, USA, 1979; pp. 34–59.

9. Spitz, K.; Moreno, J. *A Practical Guide to Groundwater and Solute Transport Modeling*; John Wiley and Sons: New York, NY, USA, 1996; pp. 120–148, 341–354, 367–379.

10. Gaiolini, M.; Colombani, N.; Busico, G.; Rama, F.; Mastrocicco, M. Impact of Boundary Conditions Dynamics on Groundwater Budget in the Campania Region (Italy). *Water* **2022**, *14*, 2462. [CrossRef]

11. Healy, R.; Scanlon, B. Groundwater recharge. In *Estimating Groundwater Recharge*; Cambridge University Press: Cambridge, UK, 2010; pp. 1–14. [CrossRef]

12. Kinzelbach, W. *Numerische Methoden zur Modellierung des Transports von Schadstoffen im Grundwasser*; R. Oldenbourg Verlag GmbH: München, Germany, 1987; p. 317.

13. Lerner, D.N.; Issar, A.S.; Simmers, I. *Groundwater Recharge: A Guide to Understanding and Estimating Natural Recharge*; IAH International Contribution to Hydrogeology: Hannover, Germany, 1990; p. 345.

14. De Vries, J.J.; Simmers, I. Groundwater recharge: An overview of processes and challenges. *Hydrogeol. J.* **2002**, *10*, 5–17. [CrossRef]

15. Scanlon, B.R.; Healy, R.W.; Cook, P.G. Choosing appropriate techniques for quantifying groundwater recharge. *Hydrogeol. J.* **2002**, *10*, 18–39. [CrossRef]

16. Lubczynski, M.W.; Gurwin, J. Integration of various data sources for transient groundwater modeling with spatio-temporally variable fluxes–Sardon study case, Spain. *J. Hydrol.* **2005**, *306*, 71–96. [CrossRef]

17. Bakker, M.; Bartholomeus, R.P.; Ferré, T.P.A. Groundwater recharge: Processes and quantification. *Hydrol. Earth Syst. Sci.* **2013**, *17*, 2653–2655. [CrossRef]

18. Gumula-Kawęcka, A.; Jaworska-Szulc, B.; Szymkiewicz, A.; Gorczewska-Langner, W.; Pruszkowska-Caceres, M.; Angulo-Jaramillo, R.; Šimůneket, J. Estimation of groundwater recharge in a shallow sandy aquifer using unsaturated zone modeling and water table fluctuation method. *J. Hydrol.* **2022**, *605*, 127283. [CrossRef]

19. Zhang, G.H.; Fei, Y.H.; Shen, J.M.; Yang, L.Z. Influence of unsaturated zone thickness on precipitation infiltration for recharge of groundwater. *J. Hydraul. Eng.* **2007**, *38*, 611–617.

20. Ounaïes, S.; Schäfer, G.; Trémolières, M. Quantification of vertical water fluxes in the vadose zone using particle-size distribution and pedology-based approaches to model soil heterogeneities. *Hydrol. Process.* **2013**, *27*, 2306–2324. [CrossRef]

21. Huo, S.; Jin, M.; Liang, X.; Lin, D. Changes of vertical groundwater recharge with increase in thickness of vadose zone simulated by one-dimensional variably saturated flow model. *J. Earth Sci.* **2014**, *25*, 1043–1050. [CrossRef]

22. Min, L.; Yanjun, S.; Hongwei, P. Estimating groundwater recharge using deep vadose zone data under typical irrigated cropland in the piedmont region of the North China Plain. *J. Hydrol.* **2015**, *527*, 305–315. [CrossRef]

23. Boughanmi, M.; Dridi, L.; Hamdi, M.; Majdoub, R.; Schäfer, G. Impact of floodwaters on vertical water fluxes in the deep vadose zone of an alluvial aquifer in a semi-arid region. *Hydrol. Sci. J.* **2018**, *63*, 136–153. [CrossRef]

24. Desaulniers, D.; Cherry, J.A.; Fritz, P. Origin, age and movement of pore water in argillaceous quaternary deposits at four sites in southwestern Ontario. *J. Hydrol.* **1981**, *150*, 231–257. [CrossRef]

25. Gehrels, J.C.; Peeters, J.E.M.; De Vries, J.J.; Dekkers, M. The mechanism of soil water movement as inferred from O¹⁸ stable isotope studies. *Hydrol. Sci.* **1998**, *43*, 579–594. [CrossRef]

26. Nolan, B.T.; Healy, R.W.; Taber, P.E.; Perkins, K.; Hill, K.J.; Wolock, D.M. Factors influencing groundwater recharge in the Eastern United States. *J. Hydrol.* **2007**, *332*, 187–205. [CrossRef]

27. Mattei, A.; Barbecot, F.; Goblet, P.; Guillon, S. Pore water isotope fingerprints to understand the spatiotemporal groundwater recharge variability in ungauged watersheds. *VZJ* **2020**, *19*, e20066. [CrossRef]

28. Cordova, J.H. Etude hydrodynamique de la Nappe Phréatique Rhénane du Secteur Saint-Louis-Bartenheim-La-Chaussée. Simulation sur Modèle Mathématique. Ph.D. Thesis, Université Louis Pasteur, Strasbourg, France, 1975; p. 126.

29. Kinzelbach, W.; Schäfer, G. Grundwassermodell S-Feuerbach, Teilbericht 2: *Numerische Untersuchungen zur Grundwassersanierung*; Universität Stuttgart: Stuttgart, Germany, 1985; p. 51.

30. McLaughlin, D.; Townley, L.R. A reassessment of the groundwater inverse problem. *Water Resour. Res.* **1996**, *32*, 1131–1162. [CrossRef]

31. R.A.S. (Regione Autonoma della Sardegna). 2005. Deliberazione della Giunta Regionale n.1/12 del 18/01/2005 con cui la Regione Sardegna ha Designato, Quale Zona Vulnerabile da Nitrati di Origine Agricola (ZVN), una Porzione del Territorio del Comune di Arborea, Cagliari. Available online: http://www.regione.sardegna.it/documenti/1_24_20050124153157.pdf (accessed on 15 February 2023).

32. Biddau, R.; Cidu, R.; Da Pelo, S.; Carletti, A.; Ghiglieri, G.; Pittalis, D. Source and fate of nitrate in contaminated groundwater systems: Assessing spatial and temporal variations by hydrogeochemistry and multiple stable isotope tools. *Sci. Total Environ.* **2019**, *647*, 1121–1136. [CrossRef] [PubMed]
33. Czekaj, J.; Jakóbczyk-Karpierz, S.; Rubin, H.; Sitek, S.; Witkowski, A.J. Identification of nitrate sources in groundwater and potential impact on drinking water reservoir (Goczałkowice reservoir, Poland). *Phys. Chem. Earth* **2016**, *94*, 35–46. [CrossRef]

34. European Community. Council directive 91/676/EEC concerning the protection of waters against pollution caused by nitrates from agricultural sources. In *Official Journal of the European Communities* L 375; Office for Official Publications of the European Communities: Luxembourg, 1991.

35. Demurtas, C.E.; Seddaiu, G.; Ledda, L.; Cappai, C.; Doro, L.; Carletti, A.; Roggero, P.P. Replacing organic with mineral N fertilization does not reduce nitrate leaching in double crop forage systems under Mediterranean conditions. *Agric. Ecosyst. Environ.* **2016**, *219*, 83–92. [CrossRef]

36. Ghiglieri, G.; Carletti, A.; Da Pelo, S.; Cocco, F.; Funedda, A.; Loi, A.; Manta, F.; Pittalis, D. Three-dimensional hydrogeological reconstruction based on geological depositional model: A case study from the coastal plain of Arborea (Sardinia, Italy). *Eng. Geol.* **2016**, *207*, 103–114. [CrossRef]

37. Pittalis, D.; Carrey, R.; Da Pelo, S.; Carletti, A.; Biddau, R.; Cidu, R.; Celico, F.; Soler, A.; Ghiglieri, G. Hydrogeological and multi-isotopic approach to define nitrate pollution and denitrification processes in a coastal aquifer (Sardinia, Italy). *Hydrogeol. J.* **2018**, *26*, 2021–2040. [CrossRef]

38. Cau, P.; Paniconi, C. Assessment of alternative land management practices using hydrological simulation and a decision support tool: Arborea agricultural region, Sardinia. *Hydrol. Earth Syst. Sci.* **2007**, *11*, 1811–1823. [CrossRef]

39. Matzeu, A. Approccio Metodologico per la Gestione Delle RIS con Particolare Riferimento Alle Aree Costiere Antropizzate. Caso di studio: Piana di Oristano. Ph.D. Thesis, University of Cagliari, Cagliari, Italy, 2018.

40. Nguyen, T.P.L.; Seddaiu, G.; Roggero, P.P. Hybrid knowledge for understanding complex agri-environmental issues: Nitrate pollution in Italy. *Int. J. Agric. Sustain.* **2014**, *12*, 164–182. [CrossRef]

41. Diersch, H.J.G. FEFLOW Software—Finite Element Subsurface Flow and Transport Simulation System, Release 5.3; Wasy GmbH Institute for Water Resources: Berlin, Germany, 2005.

42. Bear, J. *Dynamics of Fluids in Porous Media*; Elsevier: New York, NY, USA, 1972.

43. SardegnaGeoportale. Available online: https://www.sardegnageoportale.it/ (accessed on 12 July 2021).

44. Civita, M. *Idrogeologia Applicata e Ambientale*; Technical report; Casa Editrice Ambrosiana: Milano, Italy, 2005; p. 794. ISBN 8840812970.

45. Jacob, D.; Petersen, J.; Eggert, B.; Alias, A.; Christensen, O.B.; Bouwer, L.M.; Braun, A.; Colette, A.; Déqué, M.; Georgievski, G.; et al. EURO-CORDEX: New high-resolution climate change projections for European impact research. *Reg. Environ. Change* **2014**, *14*, 563–578. [CrossRef]

46. Sessini, A.M.; Ghiglieri, G.; Schäfer, G.; Lincker, M.; Cau, P.; Carletti, A. Groundwater management through Forested Infiltration Area and hydrogeological numerical modelling of Arborea coastal plain (Sardinia, Italy). In *Proceedings of the National Meeting on Hydrogeology, Napoli, Italy*, 1–3 December 2021. Available online: https://iris.polito.it/retrieve/e384c434-4f20-d4b2-e053-9f05fe0a1d67/Conference_Proceedings_Flowpath2021.pdf (accessed on 15 June 2023).

47. MENAWARA. Available online: https://www.menawara.eu/ (accessed on 24 June 2023).

48. NATMed. Available online: https://natmed-project.eu/ (accessed on 15 January 2024).

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.