Hydrogeological insights and modelling for sustainable use of a stressed carbonate aquifer in the Mediterranean area: From passive withdrawals to active management Journal of Hydrology: Regional Studies

Study area: Venafro Mts., southern-central Study focus: Via a collection of geological and hydrogeological data, a flow conceptual model of a carbonate aquifer has been coupled with a numerical model via MODFLOW code and Unsaturated Zone Flow (UZF) package in steady state and transient conditions. Simulation is further implemented with different management scenarios, for facing possible emergencies due to recharge decrease, also simulating a drastic water abstraction cut-off. New hydrological insights for the region: Carbonate fractured aquifers are a strategic water resource in the whole Mediterranean area, supplying major metropolitan areas. Despite these huge extensions, such groundwater systems are threatened by increasing drought occurrence and sig- nificant human water abstraction. A characterization of a carbonate fractured aquifer (370 km 2 ) located in central-southern Italy has been performed. Venafro Mts. Aquifer (VMA) hosts a stra- tegic resource for the Western Campania Waterworks (WCW) that supplies the populous metropolitan area of Naples, with 3.8 million inhabitants. VMA shows a slow response, with recovery time estimated at the decennial scale, testifying its limited resilience to natural and human pressures. A shift is proposed from passive management to a more comprehensive concept of smart-water monitoring, applied not only to waterworks and pipelines, but also to groundwater resources in the environment.


Introduction
Groundwater sustainability is a widespread challenge that must compete with water-well abstraction and pollution (Petts et al., 1999;Currell et al., 2012;Wu et al., 2020) as well as an increasing water demand of cities and metropolises. Urban areas of central-southern Italy are supplied by massive groundwater resources from the fractured carbonate aquifers of the Apennine range. Water springs associated with these basins have a steady flow rate and excellent hydro-chemical quality (Petitta, 2009;Mastrorillo et al., 2019), representing a strategic resource for the country. Thanks to low management costs and the vicinity of the urban areas, groundwater has facilitated the increasing water demand of Italy after the Second World War, currently 7.7 billion m 3 /year from groundwater (ISTAT database, last access 1st September 2020), favoring the 1960's economic boom, industrialization and improvement of living conditions. Economic growth and wealth increased at the same pace with the hydrogeological knowledge of these large and complex groundwater systems that encompass 17 % of Italy (Fig. 1a), with first hydrogeological works published at the end of the economic boom (Boni and Parotto, 1969). The central to southern Apennine range consists of carbonate aquifers (Boni et al., 1986) fractured and karstified, laterally limited by low-permeability boundaries. In central Italy, quantitative approaches based on a multi-year analysis of spring flow rates (Boni et al., 1986, and reference therein) identified conceptual flow models and solved hydrogeological budgets. Among these aquifers, the Simbruini-Ernici-Cairo Mts. and VMA (Venafro Mts. Aquifer) are characterized by massive groundwater resources (with a total flow rate of 27 m 3 /s) flowing into Cassino plain and Volturno valley (Fig. 1b). At the end of the 1960's, these resources were individuated to supply the rising water demands of the populous northern Campania region and the Naples metropolitan area.
Between the 1970s and 1990s, the WCW (Western Campania Waterworks) were created to exploit the groundwater resources of Cassino plain and Volturno valley, after detailed hydrogeological studies (Celico, 1983;Boni et al., 1986). However, the hydrogeological conceptual models of these groundwater basins were debated at that time and aquifer extensions were undefined.
The WCW was designed with a controversial groundwater conceptual model in support and only recently updated (Lancia et al., 2018;Saroli et al., 2019); thus, the evaluation of sustainable groundwater abstraction rates was difficult at waterworks times. Since the Fig. 1. a) Italian carbonate aquifers (in pink), a strategic resource to satisfy the national water demand; b) carbonate aquifers of central-southern Italy, red dashed line frames the VMA study area. Key to the legend: Atina gauge station (1); major springs (2); groundwater flow direction (3); WCW (Western Campania Waterworks) trace (4). c) monthly rainfall at the Atina gauge station (years 2011-2018), note the prolonged drought period during years 2012 and 2017 (c).
2000's, the groundwater resources of the Volturno Valley shortened, just several years after the WCW exploitation. Some high-elevation springs had run dry since the 2000's, and prolonged drought periods during the years 2012 and 2017 (Fig. 1c) stressed the VMA further.
However, alternative water resources are not available in the area. Increasing water demands and urbanization already put pressure on coastal and alluvial aquifers, which are damaged by pollution, sea-water intrusion and over-pumping (Petitta et al., 2013;Guyennon et al., 2017;Boumaiza et al., 2020). The central-southern Apennines carbonate aquifers, naturally safeguarded from pollution and anthropization by favorable geological and topographical conditions, are of the utmost importance and a sustainable exploit is highly desirable (Fiorillo et al., 2015).
VMA can be considered a paradigmatic case study throughout the carbonate aquifers of the central-southern Apennines, and more generally in the Mediterranean area, because of the coexistence in the same aquifer of different withdrawal waterworks, more complex than classical water supply tapping. Aquifers are frequently characterized by excessive groundwater abstraction and quality deterioration, due to high population rates, coupled with significant global changes heavily affecting groundwater resources (MED-EWI, 2007;Custodio, 2010;Romanazzi et al., 2015). Therefore, compelling management is required considering the enhanced evapotranspiration and reduced rainfall due to global warming (Giorgi and Lionello, 2008). Alarming scenarios predict a decrement of Mediterranean freshwater resources between 2-15 % for each 2 • C of warming (Cramer et al., 2018).
With the support of numerical simulations in MODFLOW with an alternative use of the Unsaturated Zone Flow (UZF) package, a modern methodology has been developed to assess the impacts of the average water abstraction performed by WCW on the spring flow rates in steady-state and transient conditions (1); in addition, it is redefined as the sustainable rate of groundwater pumping, assessing recovery time after prolonged droughts (2).

Material and methods
VMA is a minor ridge of the central Apennines with an average elevation of 700 m a.s.l. and several peaks reaching 1000 m a.s.l. (Fig. 2). This area has a Mediterranean climate with wet winters alternating to hot and dry summers. Atina gauge station ( Fig. 1c) records an average rainfall rate of 1313 mm/y between 2003 and 2018 (Regione Lazio rainfall database, last access 1st September 2020). In the last two centuries, meteorological data of southern Italy suggests a progressive rainfall depletion, especially after 1980 (Polemio and Casarano, 2004); climate has become warmer and drier, with both heavy precipitation events and long dry spells (Polemio and Casarano, 2008). Via a collection of literature works (Devoto, 1965;Celico, 1983;Saroli et al., 2014;Lancia et al., 2018;Saroli et al., 2019), geological and hydrogeological analysis allowed us to build a sound hydrogeological conceptual model, supported by a groundwater flow scheme and synthetic budgets. A conceptual model represents the basis for the groundwater numerical analysis. The numerical approach describes the ground flow system and simulates the groundwater dynamics under different conditions. Throughout the domain, a numerical model implements characteristics of the VMA and water abstraction performed by WCW. A characterization of the study site and a description of the WCW setting follow.

Geological and hydrogeological setting
The aquifer has an ideal triangular shape, bounded by tectonic elements. The study site (Fig. 2) is located at the intersection of the Roveto Valley-Atina-Caserta tectonic line on the West and the Ortona-Roccamonfina tectonic line on the East (Funicello et al., 1981). Northward, the system is interrupted by the Val Comino thrust. The stratigraphy consists of about 2000 m of carbonate deposits; basal Liassic dolostone follows Mesozoic to Cenozoic platform to basin limestone (Mostardini and Merlini, 1986). During the Apennines orogenesis, Mesozoic carbonate sequence superimposed on the Tortonian clayey sequence with the propagation of folds and thrusts elaborated by strike-slip tectonics (Centamore et al., 2007). Quaternary normal faulting fills depressions with continental to volcanic reworked material sequences (Peccerillo, 2005).
The Mesozoic to Cenozoic carbonate sequence is fractured and karstified and is schematized as a continuous aquifer bounded by less permeable Tortonian to Quaternary units. Proposed cross-sections represent the geological setting and a stratigraphic column is also included (Fig. 2). Sequences can be divided in five different technical units, according to their hydrogeological characteristics: Table 1 Elevation and average flow rate from Boni et al. (1986), Fig. 2  -DU (Dolostone Unit) is a fractured medium with good hydraulic conductivity. Tectonic activity and further weathering processes create a whitish silt corresponding with fractures and fissures. Silt hinders the groundwater flow and limits the development of a karst network (Zalaffi, 1968). Erosion carves the bedrock reaching the groundwater table. Streambed seepages are observed at Rapido River (Fig. 2). Pumping tests at Cinquina water-well ( Fig. 2) and previous numerical analysis suggest a hydraulic conductivity about equal to 10 − 5 m/s (Lancia et al., 2018). -PU (Platform Limestone Unit) outcrops at the southern sector of the domain. It consists of Mesozoic and Cenozoic limestone arranged in massive to metric layers. Under tectonic stresses, this unit has a brittle behavior (Boni et al., 1986) and fractures are distributed in accordance with the main tectonic stress. Fracture density and associated karst elements such as doline and polje drive the groundwater infiltration from the ground level to the aquifer. Pumping tests throughout the Cassino basin and previous numerical analysis suggest a high hydraulic conductivity of around 10 − 3 m/s (Lancia et al., 2018). -PBU (Platform to Basin Limestone Unit) encompasses the eastern and northern portion of the structure. This unit consists of an association of Mesozoic to Cenozoic platform, as well as scarp and basin facies (Accordi and Carbone, 1988). White and massive limestones alternate with marly limestones and chert nodules. The marl component displays a plastic behavior as the unit is subjected to tectonic stress. Plastic behavior reduces the fracturing density and slows the karst processes. However, at high elevations, karst elements such as doline and tectonic karst depressions are observed. The unit represents a regional aquifer with hydraulic conductivity estimated between 10 − 3 and 10 -4 m/s, lower than PU but higher than DU. -CU (Clayey Unit) is the regional aquitard. This unit hinders the groundwater flow thanks to its low permeability and plastic behavior under tectonic stresses. -AU (Alluvial to Volcanic Quaternary Unit) is less permeable compared to the carbonate units and consists of a multilayer aquifer (Viaroli et al., 2016). Occasionally, AU drains groundwater from the carbonate units (DU, PU, PBU) through coarser layers, if lateral contacts between the units occur.
Due to the lower hydraulic conductivity, DU has a groundwater divide role inside the VMA. The unit drains groundwater flow towards the north-western sector, whereas it feeds B1 springs (Rapido streambeds), Table 1. Nearby Valleuce village (Fig. 2), contractional tectonics with CU embedded acts as a hydraulic barrier for groundwater flow. The northern portion of Cifalco Mts. also seeps toward the DU unit and reaches the B1 springs. Conversely, the southern sector feeds S1 springs (Salauca Springs). PBU also feeds S2 springs (S. Maria d.O. Springs) and S3 springs (S. Bartolomeo Springs). S2 and S3 springs, located at higher elevations, represent the overflow levels of the VMA, as confirmed by Boni et al. (1986) recording their highly variable flow rate. B2 (Peccia streambeds) are the major springs of VMA. Thanks to carbonate superimposition at S. Pietro village, between Lungo and Cavallo Mts. (Fig. 2), B2 springs receive groundwater flows from the Venafro Mts. and additional increments from Lungo and Camino Mts. sectors. Table 1 synthetizes spring elevations and their natural average flow rate before the exploitation. For a detailed spring census, refer to Capelli et al. (2012).

The Western Campania Waterworks
Water demands of Naples and Northern Campania region (total habitants about 3.8 million) became huge after the postwar period. In addition to the historical source points, water resources of Maggiore and Matese Mts. (Fig. 1b) were exploited, causing discontent among farmers and industrialists. Thus, CASMEZ agency (Cassa per il Mezzogiorno -Development fund for the South of Italy) financed the WCW (Western Campania Waterworks) by providing 35 km of tunnels and a total length of 96 km. WCW intercepts groundwater from Simbruini-Ernici-Cairo Mts. Aquifer and VMA, via surface water intakes, well fields, and tunnel drainages. The average tapped amount is 200•10 6 m 3 /year. GW (Gari waterworks) exploited Gari Springs (Fig. 2) via shallow water intakes (Table 2). TLP (Trocchio lifting point) pumps water towards Naples and the Campania region. At Sammucro Mt. (Fig. 2), T1 waterworks consist of a tunnel with water-wells that increment the WCW flow rate. The tunnel perforates the carbonate aquifer just below the regional groundwater table and is equipped with 19 water-wells drilled inside the tunnel. Water-wells modulate the water abstraction in function of the water demand. In Venafro city, T2 tunnel drainage is equipped with lateral sluice gates that control the outflow rate. (Fig. 2). WCW exploits additional aquifers southward out of the investigated area. Sammucro and Venafro drainage systems intercept the groundwater at higher elevation with respect to natural springs, to limit operational costs. From the managerial prospective, this approach was preferred to balancing the high-fee of GW waterworks (Table 2), which requires lifting the water up to TLP (Trocchio lifting point). To overcome the hydraulic head between GW and TLP, a flow rate of 1 m 3 /s of water requires a power of 1450 kW, roughly corresponding to 0.7 million Euro per year. Waterworks characteristics of WCW are summarized in Table 2.

. Steady-state conditions
Geological and hydrogeological information as the WCW concession flow rate has been imported in a MODFLOW groundwater model and elaborated with a GMS® graphic interface. The model is solved in steady-state and transient conditions and solutions are obtained via an NWT-UPW configuration (Newton Solver-Upstream Weighting) that allows better convergences for mountainous areas (Niswonger et al., 2011). The domain consists of a single layer, discretized via 8562 square cells with a side length of 200 m. The top elevation model is imported from a STRM-DEM, resolution of 30 m. The bottom elevation is set at − 400 m a.s.l. whereas groundwater flow is essentially null (Lancia et al., 2018). The model is simplified considering the vertical lithological contacts.
A no flow boundary is set aside at the bottom as VMA is an independent carbonate aquifer ). An imposed recharge flow (RCH package) constrains the domain top and simulates the infiltration from rainfall. The model takes advantage of the Unsaturated Zone Flow (UZF) package, ordinally applied to simulate unsaturated zone flow above the water table. Neglecting the aquifer's unsaturated zone properties, the package allows simulating saturation excess and spring flow rate model as already applied for fens (Feinstein et al., 2019) and urban watersheds . The method is a valid alternative to blanketing the model with land-surface drains or other sink packages. Main springs and streambeds are represented as drains (DRN package) that allow groundwater flowing out from the domain as overpressure occurs (Fig. 3b). Initial drain conductances were estimated from field surveys (Lancia et al., 2018;Saroli et al., 2019) and further calibrated. Minor springs singularly represent less than 0.5 % of the entire water budget and are neglected for modelling scopes. The Valleluce tectonic line, where Tortonian CU are embedded into the carbonate units (PBU and DU), is represented with a barrier (HFB6 package) to locally decrease the permeability (Fig. 3b). Hydraulic characteristics of the barrier is further calibrated.
Thanks to the regional scale of the simulation, fissured to karstified bedrock is schematized as porous equivalent medium, as previously modelled in Lancia et al. (2018). Bedrock is assumed as a porous equivalent medium with an average porosity of 0.3 for the limestone and 0.2 for the dolostone. Permeability and recharge datasets are calibrated considering a single value for each unit. Calibration reduces the difference between the observed spring flow dataset and the simulated ones. Initial recharge values derive from the central Apennines budgets (Boni et al., 1986), hydraulic conductivities from pumping tests. The model uses PEST (Parameter Estimation), a nonlinear parameter estimation algorithm (Doherty, 2004) and refine the parameters. As observation dataset, calibration process considered the flow rates of the main springs observed in the 1970's and 1980's, before the water abstraction performed by WCW. For each spring, residual after calibration is lower than 0.03 m 3 /s. The impact of Sammucro and Venafro tunnel drainages on investigated groundwater flow is simulated via Water-wells (WEL1 package). Cone of depression induced by tunnel drainages is encompassed via MODPATH (Pollock, 2016). Despite the complexity of the WCW, waterworks abstraction has been schematized as sinking cells.

Transient conditions
The calibrated steady-state solution was further improved under transient conditions, between 2011 and 2018. In the simulated period, 97 stress periods of 1 month with a single time step were considered. The first step (December 2010) was set with parameters coming from the steady-state solution. Hydraulic conductivity and porosity are unit properties and only the aquifer recharge is considered as time-dependent. As direct measurement of the recharge is complex in fractured to karstified aquifers, the dataset is evaluated starting from the rainfall trend measured at the Atina gauge station (Regione Lazio database), between 2011 and 2018 (Fig. 1c). The monthly rainfall rate is processed via a convolution with a boxcar kernel considering a 3-month window. The obtained trend is further normalized according to the following formula: Where: RCH t is the recharge set at the time step t; RCH ss is the recharge calibrated in steady-state conditions; R t is the convoluted rainfall rate at the time step t; R − is the rainfall average (2011-2018) of the convoluted rainfall dataset. In the transient simulation, water withdrawals (pumping rate of WEL1 package) is considered as a steady-state. For additional simulation scenarios that consider a pumping rate cut-off, abstraction rate is homogeneously reduced from the first transient time-steps (January 2011).

Steady-state numerical results
The steady-state numerical analysis is congruent with the groundwater setting established by literature. After the calibration process, the model quantitatively reproduces the average groundwater flow observed at the major springs (Table 1). Calibrated hydraulic conductivities are in accordance with the dataset of Celico (1983) and pumping tests from unpublished reports of CASMEZ agency. The recharge dataset agrees with the theoretical value predicted by Boni et al. (1986) for central Apennines. At Valleluce divide (Fig. 3), the calibrated hydraulic conductivity of the barrier confirms its role as a hydraulic barrier of the embedded CU (Clayey Unit). In fact, S1 and B1 springs (Fig. 4a) end up being fed by different sub-basins. Table 3 reports calibrated hydraulic conductivity and recharge values.
Simulated groundwater head ranges from 210 m a.s.l. at the northern boundary to 25 m a.s.l. at B3 springs (Fig. 4a). The average hydraulic gradient is about 5.8‰, in line with the literature value for central Apennines carbonate aquifers (Celico, 1983;Boni et al., 1986). Between Cavallo and Lungo Mts. (Fig. 4a), a simulated gradient of 6.7 % is observed, whereas contact between the DU (Dolostone Unit) and PU (Platform Limestone Unit) occurs. The difference of hydraulic conductivity between the units justifies the anomaly (Celico, 1983).
The B1 springs recharge area encompasses the north-western sector of VMA; and the B2 springs the central and southern sector. In accordance with the literature review, groundwater flows from the internal sector of VMA toward B2 springs (Fig. 4b). Eastward, S2 and S3 springs recharge area covers the eastern sector of VMA. The simulated velocity is low (<0.0001 m/s) throughout the domain with peaks higher than 0.0023 m/s (Fig. 3b) close to the contacts between Dolostone Unit (DU) and Platform Limestone Unit (PU) (Fig. 2).
Numerical models elucidate the impacts of WCW on VMA in steady-state conditions (Fig. 4c). Groundwater abstraction lowers the water table throughout the domain, up to 5 m between the Montaquila and Cardito villages. T2 waterworks, with a withdrawal rate of 0.9 m 3 /s, diverges the groundwater flow originally directed at S2 springs. This spring records flow decreases of 0.95 m 3 /s, from 1.28 m 3 /s to 0.33 m 3 /s, equal to 68 % of its natural mean flow rate. At the same time T1 waterworks located in the middle of the VMA intercepts the groundwater flow directed to B2 springs. Spring flow decreases to about 0.75 m 3 /s, corresponding to 13.6 %.
According to the numerical analysis, S2 springs are affected by a lowering of 97 % of the original flow, and (Fig. 4c), as is expected from field observations and historical data, mostly dried out. The remaining flow rate (0.12 m 3 /s) is drawn from B1 and S1 springs as the water extraction triggers a diffuse water table decrease. The simulated velocity dataset slightly increases (Fig. 4d).

Transient numerical result
Transient analysis is a relevant step forward to describe the groundwater dynamics of the investigated aquifer during extended drought periods or to simulate different water abstraction rates. During the simulated period (from 2011 to 2018), the groundwater head ranges ±20 m a.s.l. from the steady-state solution. Video of the groundwater head trends under different conditions is included (Online Resource Video1.m4v).
The simulated flow rate for each spring is represented in Fig. 5. The spring flow rates following the wet and dry season cycle, with a proportional groundwater head oscillation throughout the aquifer. Minimum flow rate occurs between September and October, depending on the arrival of the first autumn rainfalls. Maximum flow rate is recorded between March and April. Compared to the bigger carbonate aquifers of the central Apennines that records maximum a flow rate between May and June (Fiorilllo et al., 2015), the maximum flow rate rate of VMA is untimely. This is due to the scarcity of snowfall, because of the limited elevation of the mountain peaks. The model also correctly simulates the 2012 and 2017 severe droughts. In this span of time, the spring flow rate shows a flattened trend, without the typical flow rate peak expected at the end of the winter in natural conditions (Fig. 5).
Introducing the water abstraction performed by WCW, B1 and B2 springs have a homogeneous drawdown during the simulated  years. According to the model, B1 and B2 springs contribute 30 % and 6 % respectively of the total groundwater abstraction performed by WCW. S1 springs do not show flow rate variations, thanks to the role of the Valleluce barrier. The remaining 64 % of the abstracted water is charged to the S2 and S3 springs. S2 springs, which quickly dry out, are penalized by the water pumping. According to the simulated dataset, groundwater flows out only during 2014 and 2015. S3 springs, also decrease to a minimum of 0.30 m 3 /s but is more stable and less affected by the seasonal trend. As S2 springs dry out, S3 springs progressively also include the recharge area of S2 springs.

Towards sustainable management of the WCW
Simulation analysis suggests VMA is severely affected by stress enhanced during drought conditions. A pumping limitation from T1 and T2 waterworks would be a possible solution to ensure an active flow of B2, S2 and S3 springs. For the same simulation period (from 2011 to 2018), further scenarios simulate the groundwater dynamics of VMA in case of a water abstraction cut-off. An initial scenario preview of T1 and T2 waterworks show a pumping rate decrease of 25 %. A second simulation introduces a drastic decrease abstraction rate of 50 %. The resulting spring flow rates are shown in Fig. 5; spring flow increase triggered by the pumping reduction are in Fig. 6. According to the model, pumping rate decrement triggers a flow rate increase of B2, of about 0.096 m3/s and 0.19 m3/s, respectively for the first and the second scenarios (Fig. 5). Located close by the T1 waterworks, S3 springs recover 0.1 m 3 /s after three months, about 0.162 m 3 /s after eight years, in the first scenario. These values rise to 0.2 m 3 /s and 0.324 m 3 /s for the second scenario (Fig. 6). Numerical analysis suggests a positive effect also for S. Maria d. O. spring, which started recovering after the 2017 drought period. After eight years of withdrawal reduction, the flow rate increase of the entire VMA corresponds to 48 % of the abstraction cutoff for the first scenario and reaches 51 % for the second one. This finding highlights the limited and slow response of VMA to management changes, due to its hydrogeological complexity. Consequently, the VMA can be treated as a fragile groundwater system, which must coexist with historical withdrawals.
New prolonged drought periods as observed in 2012 and in 2017 could further deplete the VMA groundwater resources. At the same time, the occurrence of drought periods triggered additional human needs, requiring an unavoidable full-capacity withdrawal by WCW. Consequently, the sole possible solution for limiting groundwater stress and ensuring sustainable use, is to program and perform a detailed high-frequency monitoring of natural systems. Definition of the abstraction plan based on recharge and exploitation balance and recovery strategies are an essential tool to support long-term decisions. To sum up, detailed monitoring of the VMA is compelling for a fine-tuned management plan of VMA groundwater resources In the meantime, the reduction of water lost through the water supply networks, currently 46.7 % in Campania Region (ISTAT database), is an accessible option that can provide immediate relief. Additional studies to assess the risk conditions of the surrounding carbonate aquifers are also urgent. For the studied VMA, the recommendation of an improvement of monitoring activities involves not only withdrawals, which are currently measured, but also all main springs and related water tables inside the aquifer, not forgetting rainfall and recharge rate. All this information can used with the double purpose of: i) improving the numerical model reducing its uncertainties; ii) promoting and realizing an integrated resource monitoring network which must be adopted as management tool able to guide the WCW activities.

Limits and Uncertainties of the model
The numerical model simulates head distribution and groundwater flow at the springs, under different conditions. However, numerous streamlinings are introduced. The complex stratigraphic and tectonic setting of the VMA has been strongly simplified in a single layer model, to obtain a user-friendly tool for groundwater management. Faults are implemented as hydraulic barriers only if they involve the clayey sequence characterized by low permeability. Local hydraulic conductivity changes in correspondence of faults and fractures is neglected. Recharge rate is assumed homogeneous on a yearly basis. Rainfall infiltration to VMA is preferred at dolines, tectonic-karstified depressions and polje as well as impeded by strong vegetation whereas evapotranspiration is extremely high. Methods based on the land use/cover change Shakya, 2019, 2020), merged with karst features, can better describe the recharge distribution without a perceivable variation of the total groundwater budget.
The steady-state model has been calibrated by average spring dataset flow from the 1960s and 1980s, discarding possible seasonal and long-term changes. At the same time, the water abstraction by WCW was not steady during the last 20 years, as water demand changed throughout the years. The model also does not consider the limited lateral seepage from the carbonate layers towards the coarser alluvial layers of the Cassino plain or the Volturno Valley. For the same reason, domestic water-wells, often unauthorized, are additional pressures on VMA groundwater resource. Modelling this Apennine aquifer in transient conditions appears complex with numerous intrinsic uncertainties difficult to assess. Nevertheless, the obtained results represent a significant attempt aimed at verifying the possible responses of the studied aquifer to the necessary water withdrawals, in different future scenarios.

Conclusions
Analyzed VMA suffers from a generalized groundwater depletion, due to both the historical water abstraction performed by the WCW and the prolonged dry periods more frequently observed, as happened in 2012 and 2017. Overexploiting of aquifers can damage the eco-systems and can cause economic loss to the local human activities associated with the water springs. Since the future climate scenario predicts a decrease of rainfall for the region and higher temperatures, a new management approach aimed at optimizing the groundwater abstraction from VMA is recommended.
Investigated groundwater basins, as well as probably other groundwater systems of the Mediterranean basin (Bakalowicz, 2015;Fiorillo et al., 2015) have a slow response, with the recovery time estimated at the decennial scale. Despite the limited management costs, a reduction of pumping throughout the VMA is proposed to reactivate the flow rate in some springs at risk of drying out. At the same time, a detailed monitoring of spring flow rates and groundwater heads is needed to better assess the impact of the water abstraction, both in actual conditions and in worsened scenarios. An active monitoring network of both waterworks and natural conditions will be useful to promote and integrated and sustainable management of groundwater resources. Additionally, an improvement of the entire water-supply network with active detection of piping cracks and disruptions as well as checking of illegal withdrawals would minimize the water abstraction without affecting human water demand.
Indeed, central Italy waterworks were designed and constructed before the establishment of detailed quantitative hydrogeological studies. At current times, their management must take into account increased frequency of droughts, which can trigger a diffuse wateremergency throughout the main urban areas of the country, with serious implications for the economy. Detailed studies at long-term scales are more generally suggested to assess the real sustainability of water abstraction throughout the aquifers of central Italy, following the here applied methodology. The proposed model coupled with monitoring improvement surely will represent a valid management tool to tackle current and future criticalities deriving from global changes. As added value, this approach will represent a significant shift from passive management to a more comprehensive concept of smart-water monitoring, applied not only to waterworks and pipelines, but also to groundwater resources in the environment.

Declaration of Competing Interest
The authors report no declarations of interest.