Impact of Boundary Conditions Dynamics on Groundwater Budget in the Campania Region (Italy)

Abstract

Groundwater budgets and fluxes are affected by human activities and climate change. Numerical models are cost-effective tools to investigate the different components of the hydrologic cycle. In this study, a groundwater flow model of the unconfined aquifers of the Campania region (Italy) has been developed and calibrated in Processing Modflow 11, resulting in an accurate assessment of groundwater fluxes and their trends over fifteen years (2000–2015). The model was implemented using a high-resolution grid to capture small hydrogeological features such as wells and rivers and informed by time variable datasets used as boundary conditions (i.e., river and sea levels, aquifer recharge, evapotranspiration, and discharge from adjacent systems). Good calibration and validation performances were achieved for piezometric heads (R² = 0.958). A set of scenarios was developed using constant boundary conditions (i.e., constant sea-level BC, uniform extinction depth BC), and the outputs were compared, quantitively assessing differences in groundwater fluxes. Simulations pointed out that using time series to inform boundary conditions in the model does not always result in a significant change in the computed fluxes. Overall, non-uniform extinction depth was the most influential condition, while both rivers and sea level conditions barely affected groundwater budgets. In addition, results highlighted the need for an accurate estimation of spatiotemporal variations of both recharge and evapotranspiration, due to their strong seasonal variability and their massive contribution to the hydrogeological cycle. Finally, a marked increase of evapotranspiration fluxes controlled by interannual variability of precipitation and atmospheric temperatures has been quantified over the modelled period.

1. Introduction

The Mediterranean region, where exceptional levels of biodiversity occur [1], has been identified as one of the most prominent “Hot-Spots” in actual and future climate change (CC) projections [2]. Climatic disturbance in terms of increasing temperatures [3,4], global hydrological cycle disturbance [5], and sea-level rise [6,7,8] could greatly affect both groundwater systems and groundwater-dependent ecosystems. Indeed, groundwater inputs support or compensate for surface water inputs and play a vital role in maintaining healthy ecosystems [9]. Since CC has a significant effect on the hydrological cycle [10], is mandatory for water managers to plan appropriate measures to avoid negative potential impacts on water resources [11]. Most of the studies, published so far focused on the above-ground components of the hydrological cycle, analysing both historical and projected changes, while the scientific knowledge on the subsurface components, such as recharge, groundwater levels, aquifer fluxes, and groundwater quality have been only recently emerging in literature [12,13,14]. The main tools that the research community has developed to face these issues are observation networks, which allow to record the hydraulic parameters, and mathematical/physical models which simulate their evolution in time and space. Indeed, the implementation of numerical groundwater flow models is crucial to investigating groundwater-related issues by evaluating the trends of groundwater fluxes over a simulation time and thus quantifying the ongoing CC effects on the subsurface hydrological cycle [15]. Since its release in 1988, MODFLOW has become the worldwide standard for groundwater modelling because of its flexible modular structure, complete coverage of hydrogeological processes, and free availability [16]. MODFLOW can also be integrated with a geographic information system (GIS) to provide a visual environment for groundwater resources evaluation and management [17]. Application of groundwater flow models to large domains started a long time ago, in 1978, with the Regional Aquifer System Analysis (RASA) program of the U.S. Geological Survey [18], but only a few studies, so far, aimed to numerically evaluate changes in model outputs derived from constant or variable boundary conditions (BCs).

Coelho et al. [19] developed a study focusing on the influence of different BCs on numerical simulations of an unconfined aquifer in Brazil to define the most suitable BC, without considering scenarios with constant or variable BCs. Bastani and Harter [20] developed a three-dimensional transient flow and nitrate-transport model, comparing reference results to temporally upscaled models and considering annual-averaged flow and steady-state stresses across different scenarios. However, that study only focused on the concentrations of the nitrate discharging points without quantitively addressing changes in the groundwater budget. Moreover, employing proper spatial and temporal scales for both flow and transport simulations remains one of the most important challenges in numerical modelling [21,22,23]. One of the main tasks in transient groundwater modelling is the need for large amounts of hydrogeological data for both model implementation and calibration/validation; those constraints have often limited the use of transient numerical flow models to well-characterized field sites, while transient large-scale models are still rare [24]. Despite this, transient groundwater flow modelling has been recognized as one of the hydrological research fields extremely relevant for scientific progress [25]. In recent years, open datasets on hydrology and hydrogeology have been released by many research groups and international organizations, but their widespread use in numerical models still faces problems such as the lack of universally accepted guidelines or insufficient resolution [26]. To explore the new opportunities (and limitations) offered by open datasets to represent variable BCs in large-scale groundwater modelling, the unconfined aquifers of the Campania region (Italy) were selected as representative of potentially vulnerable to CC groundwater bodies. In fact, newly released open datasets on hydrology and hydrogeology can theoretically increase the model’s capability to quantify the variations of groundwater fluxes induced by spatially and temporally variable BCs. Over the last few years, several studies have focused on developing conceptual and numerical models of the Campania region (Italy) to tackle a series of groundwater-related issues [27,28,29,30], but without investigating the impact of BCs variability on the model outputs. This study aimed to fill this gap, assessing the influence of different implementation of BCs (e.g., constant versus variable sea-level BC, uniform versus variable extinction depth BC, etc.) over the computed groundwater flow budget from 2000 to 2015. This kind of numerical analysis would also help in reducing both human and computational efforts in similar settings, by inferring what datasets are worth being simplified (e.g., averaged over time and space) and which need to be treated as spatially variable and transient over time.

2. Materials and Methods

2.1. The Study Area

The Campania region, located in the southern part of the Italian Peninsula, covers an area of about km². The region embraces a wide part of the Apennines chain, including composite geological and geomorphological units. The western slope of the Apennines is characterized by graben-like coastal depressions [31] and the plains alongside the coastal strip result from the aggradation of such structural depressions, whose bed is several thousand meters lower than the Apennine massifs [32]. The lithological and structural complexity of the region controls the hydrogeological fluxes.

The domain extends from the Tyrrhenian Sea to approximately 50 km inland, delimited in the South by the Cilento mountain range (Figure 1). The study area includes different alluvial plains (e.g., Garigliano, Volturno, Sarno, and Sele plains), characterized by debris and alluvial deposits. An important part is also occupied by limestone (e.g., M. Massico and Lattari Mts.) and volcanic units caused by the aggregation of volcanic material released during historical eruptions in the area (i.e., Vesuvio, Phlegraean Fields, and Roccamonfina volcanos are located within the model domain). In addition, the Volturno plain is characterized by a vast lowland area, subjected to land reclamation since the seventeenth century, which promoted the development of agriculture and farming along with strong coastal urbanization [33]. For toponyms used in the text refer to Figure S1.

2.2. Field Characterization and Boundary Conditions

Before the implementation of the groundwater flow model, all available datasets have been collected to define the conceptual model of the study area and the BCs (

Table 1. National and worldwide open online datasets used to set up the model.

Dataset	Parameter	Repository	URL
Soil	Grain size distribution	ISRIC	https://maps.isric.org/ accessed on 17 July 2022
Geology	Aquifer thickness	ISPRA	http://sgi2.isprambiente.it/viewersgi2 accessed on 17 July 2022
Hydrogeology	Hydraulic conductivity	ISPRA	http://sgi2.isprambiente.it/viewersgi2 accessed on 17 July 2022
CORINE Land Cover	Land use	COPERNICUS	https://land.copernicus.eu/pan-european/corine-land-cover accessed on 17 July 2022
SEALEVELEURPHYL4MY008068	Sea levels	CMEMS	https://doi.org/10.48670/moi-00141 accessed on 13 June 2022
Hydrometric station dependent	Rivers stages	CFM	http://centrofunzionale.regione.campania.it/#/pages/sensori/archivio-idrometrici accessed on 17 July 2022
MYD16A3GF	Evapotranspiration	MODIS	https://doi.org/10.5067/MODIS/MOD16A2GF.006 accessed on 10 May 2022
PCRGLOBWB 2	Recharge	UU-HYDRO	https://github.com/UU-Hydro/PCR-GLOBWB_model accessed on 12 June 2022
BIGBANG40	Evapotranspiration	ISPRA	https://groupware.sinanet.isprambiente.it/bigbang-data/library/bigbang40/grids/it_evapotraspirazione_ae accessed on 17 July 2022
BIGBANG40	Recharge	ISPRA	https://groupware.sinanet.isprambiente.it/bigbang-data/library/bigbang40/grids/it_ricarica_degli_acquiferi_gr accessed on 17 July 2022

).

The regional SRTM Digital Terrain Model with a spatial resolution of 20 × 20 m cells was used to retrieve topographic information and represent morphology. The data files contain the elevation in a digital format, averaged on the model grid. Vegetation, buildings, and other cultural features were removed digitally, resulting in a proper terrain model.

The hydraulic conductivity of lithological units was obtained from permeability tests, by mining data from 91 borehole reports available on the online database of Instituto Superiore per la Protezione e Ricerca Ambientale [34]. All the data about pumping well rates, effective radius, screen length, and groundwater drawdown were retrieved from the same reports. The depth of the unconfined aquifer for each borehole was collected along with locations (latitude, longitude, and UTM coordinates). Other information on the aquifer’s depth has been retrieved from previous studies carried out in the Campania region, such as logs on the Volturno river plain [35] or aquifer information in Cuma [36,37], Bagnoli, Giuliano, and San Valentino Torio.

Measurments on sea levels were downloaded from the open dataset Servizio Mareografico Nazionale (SMN). SMN exploits a network of 36 measuring stations uniformly distributed over the national territory and located mainly inside the port facilities. The hydrometric level is monitored since 2010 using a new microwave (radar) level meter with millimetric precision, coupled with a second level meter with a floater. Each level is referred to the absolute altimetric network of Istituto Geografico Militare (IGM). To gain a better understanding of the sea level behaviour across the entire region, Napoli, Salerno, Gaeta, and Palinuro stations were considered. However, since SMN open dataset does not cover the entire period of simulation, the datasets of the Copernicus Marine Environment Monitoring Service (CMEMS) were also exploited for sea-level trends at the same locations. CMEMS altimetric satellite data were elaborated together with the SMN ones allowing to cover the entire period of study.

The Moderate Resolution Imaging Spectroradiometer (MODIS) products have been deployed to obtain the trends of evapotranspiration in the study area for the period 2000–2015. Data were downloaded using the AppEEARS [38] interface and the MYD16A3GF model [39] to obtain spatially interpolated values. The algorithm used in MODIS is based on the Penman-Monteith equation, which includes daily meteorological data along with remotely sensed variables, such as vegetation dynamics, albedo, and land cover [40]. Land use and soil texture maps (Table 1) were overlapped via GIS to derive a two-dimensional matrix representing the evapotranspiration extinction depth [41].

Time series of groundwater recharge has been extracted from the outputs of PCR-GLOBWB v2.0 (Utrecht University, Utrecht, The Netherland) [42], a grid-based global hydrology and water resources model. Global hydrological models (GHMs) are tools able to quantify the global terrestrial water cycle, simulating the distributed hydrological response to weather and climate variations at high resolution.

Hydrometric levels were downloaded from the Centro Funzionale Multirischi (CFM), a Civil Protection Agency of the Campania region, to represent the variation of the river stage over the simulation period. CFM provides a list of all hydrometric stations of the region recording data with a daily frequency. The closest stations to the model’s boundaries were considered representative of the baseflow entering the domain from the fractured aquifers. Hydraulic head fluctuation in the fractured aquifers was modified according to the average baseflow variations over the simulation time.

2.3. Numerical Modelling

To estimate the groundwater flow budget at a regional scale, assessing the influence of variable BCs, all the data were elaborated in a GIS environment and processed via Processing Modflow 11. The flow simulation was carried out using the MODFLOW-2005 code [43]. The equation that stands at the base of MODFLOW-2005 is the partial-differential equation that describes the three-dimensional movement of groundwater of constant density through a porous matrix:

$$\frac{\partial }{\partial \mathrm{x}}\left({\mathrm{K}}_{\mathrm{xx}}\frac{\partial \mathrm{h}}{\partial \mathrm{x}}\right)+\frac{\partial }{\partial \mathrm{y}}\left({\mathrm{K}}_{\mathrm{yy}}\frac{\partial \mathrm{h}}{\partial \mathrm{y}}\right)+\frac{\partial }{\partial \mathrm{z}}\left({\mathrm{K}}_{\mathrm{zz}}\frac{\partial \mathrm{h}}{\partial \mathrm{z}}\right)+\mathrm{W}={\mathrm{S}}_{\mathrm{s}}\frac{\partial \mathrm{h}}{\partial \mathrm{t}}$$

(1)

where K_xx, K_yy and K_zz are values of hydraulic conductivity along x, y and z coordinate axes (LT⁻¹); h is the potentiometric head (L); W is a volumetric flux per unit volume representing sources and/or sinks of water (T⁻¹); S_s is the specific storage of the porous material (L⁻¹); and t is time (T). The standard version of MODFLOW-2005 treating the fractured units as equivalent porous media and not complex ones such as the Conduit Flow Package [44] was used given the lack of information on the location and dimensions of such zones. Thus, for the principle of parsimony in numerical modelling, a simpler representation was here adopted [45].

The model domain consisted of 656 rows, 335 columns, and 3 layers with variable width and thickness. The top of the grid represents the local ground surface, ranging from −3.49 m to 1642 m above sea level (a.s.l.). The bottom of the grid represents the bottom of the unconfined aquifers and ranges from −83 m below sea level (b.s.l.) to 144 m a.s.l.

Initially, a steady-state simulation has been run and the model has been calibrated via PEST [46] against an observed piezometric map of the Campania region for the year 2000 [47]. Given that the piezometric map was available only as a vector map, it was sampled along the contour lines to create 500 observation points that were compared versus calculated ones with an estimated error of ±1.0 m. Tyrrhenian Sea levels have been initially set to 0 m as constant head cells. Values of hydraulic conductivity, corresponding to boreholes in the same cluster, have been averaged to get an average value for each geological formation in the study area. In this way, the model domain was subdivided into polygons whose extension corresponds to the extension of the different geological formations reported in the geological map. Figure 2 highlights the heterogeneity of the Campania region, where the horizontal hydraulic conductivity values range from 1 × 10⁻⁶ m/s in the coastal part of the Volturno plain, where clay layers are present, to 6.6 × 10⁻⁴ m/s in the Vesuvio area, where highly fractured volcanic rocks are present. The value of the vertical hydraulic conductivity was set equal to the value of the horizontal hydraulic conductivity assigned to that cell in fractured units, while it was set as a tenth of the horizontal hydraulic conductivity in sedimentary units, to account for internal layering. It must be stressed, however, that local scale heterogeneities such as intensively faulted zones, small lenses, etc. have not been explicitly incorporated in this large-scale model.

For the steady-state simulation, evapotranspiration rates for the year 2000 were averaged from the MODIS dataset and interpolated on the numerical grid; while the recharge rates were averaged from the PCR-GLOBWB v2.0 for the same period and interpolated on the numerical grid. The Processing Modflow polyline input method in the River Package was employed to simulate, from North to South, Garigliano, Savone, Agnena, Volturno, Regi Lagni, Sarno, Picentino, Tusciano, Sele, and Calore rivers, drawing their actual pathways and defining river parameters. The riverbed conductance per unit length C_rivl (LT⁻¹), the river stage h_riv (L), and the elevation of the riverbed bottom B_riv (L) have been assigned to the most upstream and downstream vertices of each river polyline and then applied to the model cells employing a linear interpolation. River stages were assumed close to the assigned starting head value and kept constant for the entire steady-state simulation. The flow rate Q_riv (L³T⁻¹) between each river reach and the model cell is calculated by the following equations:

$${\mathrm{Q}}_{\mathrm{riv}}=\left\{\begin{array}{c}{\mathrm{C}}_{\mathrm{riv}}·\left({\mathrm{h}}_{\mathrm{riv}}-\mathrm{h}\right), \mathrm{i}\mathrm{f} \mathrm{h}>{\mathrm{B}}_{\mathrm{riv}}\\ {\mathrm{C}}_{\mathrm{riv}}·({\mathrm{h}}_{\mathrm{riv}}-{\mathrm{B}}_{\mathrm{riv}}), \mathrm{i}\mathrm{f} \mathrm{h}<{\mathrm{B}}_{\mathrm{riv}} \end{array}\right.$$

(2)

where h (L) is the calculated hydraulic head in the cell containing the reach.

The wells analysed from the online ISPRA database were entered into the model through the Well Package. Industrial wells were set with continuous pumping rates, while the agricultural wells were set with discontinuous pumping rates, active only during the spring and summer seasons. Since the specified wells penetrate more than one model layer, the total pumping rate has not been distributed over the penetrated layers, but it has been assigned to the lowest cell of the well, belonging to the third model layer. The Drain Package was used to simulate the reclamation channels in the Marina di Licola area and also to simulate the effect of the underground metro lines in Naples, the mean height of each feature a.s.l. was derived from the topographical online database, while the cell conductance C_d (L²T⁻¹) was set to a very high value (0.01) to ensure optimal groundwater withdrawal from the aquifer. The location of the depth of the different metro lines was sourced from Russo et al. [48]. The discharge rate Q_d (L³T⁻¹) from the model cell to the drain is calculated by:

$${\mathrm{Q}}_{\mathrm{d}}=\left\{\begin{array}{c}{\mathrm{C}}_{\mathrm{d}}·\left(\mathrm{h}-{\mathrm{h}}_{\mathrm{d}}\right), \mathrm{i}\mathrm{f} \mathrm{h}>{\mathrm{h}}_{\mathrm{d}}\\ 0, \mathrm{i}\mathrm{f} \mathrm{h}\le {\mathrm{h}}_{\mathrm{d}} \end{array}\right.$$

(3)

The horizontal hydraulic conductivity was estimated by inverse modelling using the numerical code PEST, assigning a parameter number to it. PEST searches optimum parameter values for which the sum of squared deviations between model-calculated and observed values from the piezometric map of the region is reduced to a minimum. Such comparison has been performed using a scatter diagram of the simulated and observed data (Figure 3).

To investigate possible changes in the groundwater system associated with temporal variations of evapotranspiration and recharge patterns, a transient simulation was carried out. Sixty-four stress periods of seasonal length (3 months), subdivided into three time steps of 1 month each, were used. As recommended by Harbaugh [43], the starting hydraulic heads for the transient simulation have been set equal to the head values calculated at the end of the steady-state simulation.

Table 2. S_s and S_y depending on the geologic material (upper panel), and C_riv gained after the calibration procedure (lower panel).

Material	Ss (m−1)	Sy (-)
Clayey and calcareous clayey units	0.001	0.2
Debris, alluvial, fluviolacusters and fluvioglacial deposits	1 × 10−4	0.33
Lavas, ignibrites, pyroclastics	1 × 10−5	0.21
Organogenic, neritic, and platform limestone	1 × 10−5	0.14
River	Criv (m2/s)
Garigliano	from 0.021 to 0.235
Savone	from 0.017 to 0.058
Agnena	from 0.001 to 0.033
Volturno	from 0.009 to 0.044
Regi Lagni	from 0.007 to 0.072
Sarno	from 0.007 to 0.040
Picentino	from 0.009 to 0.011
Tusciano	from 0.011 to 0.019
Sele	from 0.014 to 0.116
Calore	from 0.060 to 0.086

shows the specific storage (S_s) and specific yield (S_y) values, depending on the geologic material, together with the riverbed conductance (C_riv) gained after the calibration procedure.

To consider the variability of the evapotranspiration over time, the MODIS data have been spatialized using the kriging technique and inserted into the model with seasonal frequency. The values of recharge rate from the global numerical hydrologic model published by de Graaf et al. [49] have been processed in python to obtain spatially interpolated values for the period 2000–2015, which have been assigned to the model with seasonal frequency. The calibrated model was implemented with the before-mentioned transient datasets accounting for the seasonal changes in hydrometric levels, Tyrrhenian Sea levels, and fractured aquifer seasonal head fluctuations as variable BCs. The Evapotranspiration Package was also integrated with spatially variable extinction depths dependent on the land use and soil texture maps.

The validation of the transient state model was performed by exploiting a piezometric map of the Campania region for the year 2013, which was derived from a secondary elaboration of the global groundwater depth dataset [50] to match modelled and measured hydraulic heads. Besides, selected monitoring wells (Figure 1) with piezometric trend data were compared versus model results.

Both the model calibration and validation were assessed by calculating and comparing the value of the coefficient of determination (R²) and the mean absolute error (MAE).

The Processing Modflow ZoneBudget routine [43] can produce water budgets for user-defined zones within the model domain. Zones are defined by assigning zone numbers to model cells, allowing for to calculation of the volumetric water budget in a different part of the domain, accounting for evapotranspiration, fractured aquifers, rivers, and the Tyrrhenian Sea component separately. Each trend calculated by the model, considering all the transient input datasets as variable BCs (transient model), was compared with the respective trend calculated by an upscaled model considering a uniform extinction depth, a constant head at the fractured aquifers, and constant sea and river levels (basic model). Moreover, to investigate which datasets are worth being averaged, a set of scenarios have been developed applying to the transient model each of the averaged BCs individually:

• constant head at the fractured aquifers BC (scenario S1);
• uniform extinction depth BC (scenario S2);
• constant Tyrrhenian Sea level BC (scenario S3);
• constant river levels BC (scenario S4).

In addition, to numerically quantify the impact of the averaged conditions, the 15-year average flux of each groundwater budget component has been estimated and differences between basic, transient, S1, S2, S3, and S4 simulations were graphically compared.

Finally, since the pivotal contribution to the groundwater budget of both recharge and evapotranspiration, an additional recharge and evapotranspiration rates dataset has been investigated. Alternatively, to the PCR-GLOBWB v.2.0 published by Sutanudjaja et al. and the MODIS products, the BIGBANG (Bilancio Idrogeologico GIS BAsed a scala Nazionale su Griglia regolare) [51] distributed model was employed, and results compared.

3. Results and Discussion

3.1. Flow Model Performance

Given the extremely challenging regional scale of the domain (3672 km²), reasonably good calibration and validation were achieved for hydraulic heads, resulting in an R² of 0.958 (Figure 3). MAE values of 3.56 m for calibration and 2.87 m for validation, were judged acceptable in the wide range of the simulated piezometric levels (>100 m). Observations of groundwater levels from selected monitoring wells were used to verify the ability of the model to describe water table fluctuations in the unconfined aquifers (Figure 3). Unfortunately, a widespread monitoring network of wells is not available in the region, and this is a limitation for the model performance. Indeed, monitoring groundwater levels at regional scale may enable a better assessment of groundwater depletion, as well as inferring changes in flow directions over time, as recently pointed out by Brakkee et al. [52]. Moreover, piezometric data are missing in some karst and fractured areas of the domain (Massico mountain, Vesuvio volcano, and Monti Lattari) consequently model reliability is lower on those areas than in alluvial zones. By means of the Processing Modflow ZoneBudget routine [43] a defined zone number was assigned to those areas and their importance in relation to the overall model domain was assessed in terms of water budget component percentages. As expected, groundwater recharge played a very important role in karst environments, representing 12.7% of the overall domain, while evapotranspiration was just 1.5% with respect to the total, given that in mountain areas the groundwater depth from the surface is higher than the extinction depth. The storage capacity turns out to be consistent, whereas all the other budget components resulted to be negligible. The general high permeability grade due to fracturing leads to a high groundwater recharge rate [53] that varies between 48% and 78% of the mean annual effective precipitation [54,55], depending on the relative abundance of limestone and dolomites lithologies.

The validated transient model allowed analysing and quantifying relevant stresses in different portions of the domain at different times. Figure 4 shows an example of both the three-dimensional and planar view of the calculated hydraulic heads map at the end of the simulation. It pointed out regions of the domain characterized by high hydraulic gradient, especially the northern part in correspondence of the Roccamonfina Volcano, which acted as a groundwater mounding area, as well as the Vesuvio, Phlegraean Fields, and Monti Lattari. Hydraulic heads resulted to be very low in the Volturno river plain, following the topographic gradient. Near the coast, the water level was very shallow, due to reclamation channels located in the Marina di Licola area (lower Volturno plain).

In addition, the model produced accurate estimates of the groundwater budget with a mean error of less than 0.3% among in and out components. This allowed focusing on the single water budget components, whose fluxes are difficult to estimate without a physically based model even though they are crucial for better management of the hydrogeological basins.

3.2. Water Fluxes Trends and Comparative Analysis

The fluxes trends have been estimated accounting for the inflow from the fractured aquifers, the evapotranspiration, the river-groundwater interactions, and the Tyrrhenian Sea component, separately. To highlight the influence of the time variable datasets, the basic model calculated trends (red line), the transient model calculated trends (blue line), and the scenarios (S1, S2, etc.) calculated trends (dashed green line) associated with the budget component have been plotted together.

Figure 5a displays the inflow component from the fractured aquifers. The transient model calculated trend shows a quite higher flux with a wider seasonal variation throughout the simulation time, with respect to the basic and S1 model.

The evapotranspiration component represents the most important outflow contribution. Figure 5b clearly shows its seasonal variability with increasing fluxes over the warm seasons, reaching values up to 20 m³/s. Moreover, since the extinction depth affects the evapotranspiration fluxes, the S2 trend perfectly reflects the basic model, underlying the importance of considering a non-uniform extinction depth over the domain to account for different land cover and different soil types that directly drive the extinction depth [41].

The seawater intrusion component (Figure 5c) is characterized by a very low flux and can be considered negligible with respect to submarine groundwater discharge. However, importing a transient sea-level dataset leads the seawater intrusion flux of the transient model to display positive peaks due to the seasonal oscillation of the Tyrrhenian Sea. This could be further investigated at the local scale by performing a variable-density flow and transport simulation. There are no evident changes between the basic, transient, and S3 scenarios in submarine groundwater discharge, which shows a mean flux over the simulation time between 5 and 6 m³/s (Figure 5d).

Both river inflow and outflow (Figure 5e,f), which represent the losing and gaining flow rate from and to the river, show no significant changes between basic, transient, and S4 simulations. Therefore, considering the effort to import all required data in the models, a constant river level dataset averaged over the simulation time should be considered for large-scale simulations.

To numerically quantify the impact of time variable input datasets on groundwater budget, the 15-year average flux of each component has been estimated and differences between basic, transient, S1, S2, S3, and S4 simulations were represented in Figure 6. A difference of ~1.67 m³/s was found in the mean evapotranspiration between the basic and the transient simulations. Almost the same difference (~1.70 m³/s) was found between S2 and the transient model. This confirmed that extinction depth is the most influencing parameter, requiring a precise spatiotemporal estimation to obtain accurate predictions of the groundwater budget.

The difference in the mean inflow from the fractured aquifers component (~0.65 m³/s) between the basic and transient models did not seem to be entirely due to the associated transient inflow parameters, but also accounted for the variable river level dataset.

Even if the difference in the averaged seawater intrusion component was relatively high (almost doubled from the basic to transient model), it should be considered negligible in absolute terms (0.02 m³/s). The mean submarine groundwater discharge showed no appreciable changes between the scenarios, both in relative and absolute terms. Effectively, that component seems to be more affected by the non-uniform extinction depth dataset (0.16 m³/s difference) than by the transient Tyrrhenian Sea levels (0.06 m³/s difference), suggesting that sea level could be kept constant over the simulation time.

The same analysis can be conducted considering the mean groundwater flux towards the river (river outflow). The difference of 0.98 m³/s between the basic and the transient simulations was mostly due to the non-uniform extinction depth rather than to the variable river level dataset: the value rose to 1.02 m³/s between the S2 and the transient model, while it dropped to 0.05 m³/s between the S4 and the transient model. Considering also that the mean river flux towards the unconfined aquifers (river inflow) slightly differs among the various scenarios, accounting for a variable river level dataset does not seem to be crucial.

The comparative analysis pointed out the fundamental role of the extinction depth in shifting the allocation of groundwater fluxes among different compartments, altering the contribution of each component with respect to the total budget.

Moreover, since the evaluation of each groundwater contribution underlines the pivotal role of both the recharge (with a mean flux of 18.1 m³/s) and evapotranspiration components, the additional BIGBANG dataset has been investigated. The derived recharge and evapotranspiration trends have been calculated and compared with the transient model results (Figure 7).

The recharge flow rates strongly differ among the two models, while the evapotranspiration fluxes were more similar. The linear temporal trends show an increase in the recharge flux over time for both the PCR-GLOBWB v2.0 and the BIGBANG products, reflecting the positive trends of the hydraulic heads observed in the monitoring wells P1, P2, P3, and P4 (Figure 3). The process of groundwater recharge is driven by a decadal variability of the North Atlantic Oscillation and its robust correlation with the rainfall patterns [56]. Both the calculated evapotranspiration fluxes, controlled by interannual variability of precipitation and air temperature, showed increasing trends with a small positive shift (~0.10 m³/s per year), highlighting that CC impact on groundwater resources is confirmed by the two independent databased here employed. Furthermore, because of the overestimation of the recharge flux due to topographic and geological heterogeneities of the Campania region, running the transient model using the BIGBANG products slightly decrease the performance parameters, with an R² = 0.948 and an MAE = 3.42 m, producing flooded areas during the peak recharge seasons that were not physically based.

Finally, it must be stressed that all other scenarios showed very little variation of the model performance indicators with an R² and an MAE that never exceeded 0.947 and 2.91 m, respectively. This is due to the low sensitivity of piezometric heads in capturing the changes induced by evapotranspiration at the basin scale [57], while to solve this complex issue, monitoring wells with continuous heads dataloggers should be employed coupled with baseflow measurements [58].

4. Conclusions

Open online datasets, along with historical observations and literature review were used to untangle the complexity of regional hydrogeological settings, supporting the implementation of a robust numerical model to describe groundwater fluxes for the unconfined aquifers of the Campania region. The quantification of the total flux of each component in the budget highlighted the pivotal role of recharge and evapotranspiration rates, with a mean flux of 18.1 and 12.3 m³/s, respectively. Their variable patterns throughout the simulation period affected the whole hydrogeological budget, pointing out the need for their accurate estimation over time, to account for their high seasonal variability and climatic trends. Recharge can reach a value of above 60 m³/s during wet periods, while evapotranspiration variability ranges from 5 to maximum values up to 20 m³/s over the warm seasons. Time variable input datasets had very different effects on the components of the groundwater budget, controlling water fluxes between different compartments. Simulations pointed out that both fractured aquifers head fluctuations and non-uniform extinction depth datasets used as BCs should not be averaged. Considering transient data of extinction depth had a massive effect on the budget results, not only influencing evapotranspiration rates directly but also affecting all the other components of the groundwater budget, indirectly. Conversely, the inclusion of other datasets, namely variable river and sea levels, barely impacted the groundwater budget, not making it worth the extra effort. Therefore, the inclusion of transient BCs in a numerical model should be based on expert judgment, balancing accuracy and complexity and considering the human and computational efforts required. Finally, the results have underlined that continuous groundwater level and baseflow estimates are crucial to minimize model uncertainties and increase the reliability of groundwater budgets at the basin scale. Therefore, such measurements need to be included in an ideal groundwater monitoring network at a regional scale.