Next Article in Journal
Evolution of Sustainable Water Resource Utilization in Hunan Province, China
Next Article in Special Issue
Simulation of Heat Flow in a Synthetic Watershed: Lags and Dampening across Multiple Pathways under a Climate-Forcing Scenario
Previous Article in Journal
Evapotranspiration Measurements and Modeling
Previous Article in Special Issue
Improving Results of Existing Groundwater Numerical Models Using Machine Learning Techniques: A Review
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Water
Volume 14
Issue 16
10.3390/w14162475
water-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Stepwise Modelling Approach to Identifying Structural Features That Control Groundwater Flow in a Folded Carbonate Aquifer System

by
Elisabetta Preziosi
1,
Nicolas Guyennon
1,
Anna Bruna Petrangeli
1,
Emanuele Romano
1 and
Cristina Di Salvo
2,*
1
CNR-Water Research Institute, Area della Ricerca di Roma 1, Via Salaria Km 29,300-C.P. 10, Monterotondo Stazione, 00015 Rome, Italy
2
CNR-Institute of Environmental Geology and Geoengineering, Area della Ricerca di Roma 1, Via Salaria Km 29,300-C.P. 10, Monterotondo Stazione, 00015 Rome, Italy
*
Author to whom correspondence should be addressed.
Water 2022, 14(16), 2475; https://doi.org/10.3390/w14162475
Submission received: 5 July 2022 / Revised: 8 August 2022 / Accepted: 9 August 2022 / Published: 11 August 2022
(This article belongs to the Special Issue Groundwater Hydrological Model Simulation)
Browse Figures
Review Reports Versions Notes

Abstract

:
This paper concerns a stepwise modelling procedure for groundwater flow simulation in a folded and faulted, multilayer carbonate aquifer, which constitutes a source of good quality water for human consumption in the Apennine Range in Central Italy. A perennial river acts as the main natural drain for groundwater while sustaining valuable water-related ecosystems. The spatial distribution of recharge was estimated using the Thornthwaite–Mather method on 60 years of climate data. The system was conceptualized as three main aquifers separated by two locally discontinuous aquitards. Three numerical models were implemented by gradually adding complexity to the model grid: single layer (2D), three layers (quasi-3D) and five layers (fully 3D), using an equivalent porous medium approach, in order to find the best solution with a parsimonious model setting. To overcome dry-cell problems in the fully 3D model, the Newton–Raphson formulation for MODFLOW-2005 was invoked. The calibration results show that a fully 3D model was required to match the observed distribution of aquifer outflow to the river baseflow. The numerical model demonstrated the major impact of folded and faulted geological structures on controlling the flow dynamics in terms of flow direction, water heads and the spatial distribution of the outflows to the river and springs.

1. Introduction

Carbonate aquifers are important groundwater resources worldwide due to their high permeability and rapid groundwater velocities. Ford and Williams [1] estimate that 20% of the world population largely depends on groundwater from carbonate aquifers. These are often defined as karst aquifers due to their “self-organized, high permeability channel networks formed by positive feedback between dissolution and flow” [2]. An important feature of the carbonate aquifers, especially when diffuse flow prevails, is their capability to store large quantities of groundwater during humid periods and gradually release them during dry periods. Hence, they are fundamental to sustain both human uses and groundwater-related ecosystems in many parts of the world. Water quality in carbonate aquifers is often excellent; hence, they are regarded as strategic both for human consumption as well as to sustain environmental uses. The high permeability often results in thick unsaturated zones, so exploitation of groundwater is often from low-elevation springs, especially in mountainous areas [3]. Pumping wells located near the springs are used to overcome spring discharge shortage in dry seasons [4].
These resources are very often exploited to supply large urban areas, and the evaluation of the possible negative effects of climate changes on their discharge is challenging. The decrease in annual precipitation and increase in temperature and evapotranspiration due to climate change in the Mediterranean area [5,6,7,8,9] is leading to a general decline both in surface water discharge [10] and in groundwater recharge [11]. Preziosi and Romano (2013) [12] highlighted that the spring discharge of large carbonate aquifers in Central Italy has decreased by 6% to 42% in the period 1938–2011, the highest percentages estimated for the smallest springs [13].
Recent advances in groundwater flow and level prediction include the application of diverse methods including wavelet analysis, Gaussian process regression [14], machine learning modes [15] and integrated numerical modeling [16].
The development of numerical models is considered a fundamental step for the adoption of water management plans aiming to preserve groundwater resources and the related ecosystems [17]. Many authors have focused on the numerical modelling of karst systems to assess the risk of spring discharge shortage due to climate change [18,19] or to evaluate the effects of withdrawals [4]. However, numerical modeling of carbonate aquifers in folded and faulted terrains is a challenge due to the complexity of the hydrogeological systems, and excessive simplification may lead to an unsatisfactory predictive capability of the model. Carbonate aquifers are often characterized by highly conductive conduit flow paths embedded in a less conductive fissured and fractured matrix [20]. In spite of these strong permeability contrasts, the equivalent porous medium approach (EPM) can be applied to karstified aquifers with some limitations; specifically, it may not be suitable for deterministically modeling flow along faults, and it often fails to predict flow direction and velocity. However, it can correctly approximate flow and spring discharge at the regional scale [21]. Further, Abusaada and Sauter (2013) [22] affirm that EPM models can simulate flow in karst aquifers as long as the simulated saturated volume is large enough to average out the local influence of karst conduits. The significant influence of the geological structure (especially folding and lithology) and the karst system on the location of the springs and their flow regime has been addressed by [23,24]. Structural folds may divert groundwater flow from the general hydraulic gradient. The presence of marl layers may sustain perched sub-aquifers above the regional aquifer, and karstification may locally increase the hydraulic conductivity by several orders of magnitude. Moreover, the complex geometry of model layers can result in different thicknesses of saturated portions, which can imply the drying and rewetting of cells during model iterations, leading to numerical instabilities, preventing convergence and increasing numerical error [25]. Including all these characteristics in the conceptual and numerical model requires the adequate definition of layer thickness, dipping and hydraulic properties. Recent advances in groundwater numerical modeling include the development of solvers which facilitate achieving convergence and/or reducing computational errors due to model nonlinearities, as well as packages tailored for solving specific problems [26]. This allows the development of fully 3D numerical models which are able to reproduce complex settings with stable solutions. However, by increasing model structure complexity, the number of input parameters increases as well as the related uncertainty. A model with too many parameters is susceptible to over-fit the data [27,28,29]. The higher the complexity, the more accurate the calibration procedure should be, requiring an adequate number of calibration targets. Simulating complexity not supported by the data can be useless and misleading.
The aim of this research was to develop and test a modelling procedure for the simulation of groundwater flow in a complex karst, folded, multilayer aquifer using the EPM approach. In this framework, three steady-state numerical models of a carbonate aquifer in Central Italy (Monte Coscerno) were developed with increasing complexity as warranted by the inability of the simpler model to adequately reproduce observations [30]. A stepwise procedure was developed for testing the ability of the models to reproduce observations. Seeking parsimony, we compare the results from a simple one-layer model and more complex quasi-3D and fully 3D models with a different number of layers and spatial variability of parameters.

2. Conceptual Model

2.1. Geological and Hydrogeological Setting

The case study aquifer (Monte Coscerno, 230 km2) is located in the Apennine Range in Central Italy (Figure 1). It can be described as a box-fold anticline with a nearly meridian axis and a vertical-to-overturned forelimb [31] belonging to the frontal thrust ramp of a regionally curved thrust known as the Mount Coscerno–Rivodutri thrust [32,33]. The above mentioned thrust on the east and the Valnerina thrust [34] on the west, with an approximately meridian direction, bound nearly all the aquifer and act as a no-flow boundary (black solid line in Figure 1). The general groundwater flow is mainly in a south to north direction, parallel to the prevailing tectonic lines [35,36,37]. The areas where the recharge is most effective are the plateaus of Monte Coscerno, Monte Aspra and other summits that exceed 1000 m a.s.l., reaching 1685 m a.s.l. at Monte Coscerno (Figure 2). A sequence of Meso-Cenozoic calcareous formations interbedded with marl layers results in a multilayer aquifer system, with three main sub-aquifers (from bottom to top: Calcare Massiccio-Corniola, Maiolica, Scaglia limestone units) separated by two aquitards [36]. The aquitards are locally discontinuous due to depositional, erosional and tectonic effects, favoring vertical leakage between the three aquifers [37,38]. The base of the aquifer is represented by the top of the Triassic dolomitic evaporitic complex (“Marne a Rhaetavicula contorta” marlstones and “Anidriti di Burano” anhydrites, [39,40], mostly dipping westward. Groundwater is expected to flow according to the bedding attitudes in the direction of the steepest structural descent. This may result in large unsaturated portions in the eastern part of the aquifer (Figure 3); moving from east to west, the sub-aquifers 2 and 3 become confined, feeding the Scaglia sub-aquifer through vertical leakance upward. The sub-aquifers 1 and 2 do not extend through all the model area due to erosional processes which affected the anticline eastside where the more ancient formations outcrop (Figure 1 and Figure 3). For this reason, rainfall infiltrates through the uppermost outcropping aquifer and flows westward and northward according to the dip of the layers.
Hourly hydrometric data and periodic discharge measurements of the Nera River in two stream gauging stations (Vallo di Nera and Torre Orsina, Figure 1) have been made available by the Umbria Region since 2006. The monthly discharge of the Lupa spring is available for the period 1985–1997. Moreover, the Lupa spring daily discharge (since 1998) and piezometric heads measured in the Scheggino well (since 2001) are available online from the local Regional Environmental Agency [41]. Head data are very scarce, except for the already mentioned Scheggino well and the very recently installed Renari di Capriglia well (Figure 1). The Nera River is incised into the carbonate aquifer, increasing its discharge from north to south. Several discharge measurements performed in the years 1991–1993 at the 6 gauging stations in Figure 1 [42] indicate a conspicuous discharge increment, nearly constant throughout the year, revealing gaining stream conditions between Vallo di Nera and Umbriano (reaches R1 to R4, Figure 1). The river baseflow was estimated between 3.2 and 3.4 m3/s in the period 1991–1993 by Boni and Preziosi (1993) [35]. In addition, the aquifer feeds several point springs (Scheggino, Lupa and Pacce) and the Precetto stream (Figure 1, Table 1).
The aquifer discharge to the Nera River was estimated as the discharge increment between two gauges (Vallo di Nera and Torre Orsina gauges, Figure 1) using spot measurements in the years 1997–2012 provided by the Regional Environmental Agency [41]. The average of the 83 spot measurements is 3.28 m3/s, ranging from 0.9 to 7.46 m3/s (Figure 4). The total aquifer discharge was estimated at about 3.4–3.6 m3/s, with an extremely regular seasonal regimen uncommon in karst areas [42]. However, there is no evidence in the area of developed karst conduits despite the presence of likely fractured limestones. Consequently, Monte Coscerno can be classified as a “diffuse flow aquifer” [44] i.e., a carbonate aquifer system with dispersive circulation due to a micro-fractured interconnected network with extremely reduced or even inexistent karstification without preferential drainage paths [45]. There is no concentration of flow towards localized springs, with the exception of the Lupa, Pacce and Scheggino springs. The Nera River flows parallel to the Valnerina thrust, as shown in Figure 1, that represents the no-flow boundary of the aquifer on the eastern side at the lowest topographic elevation, explaining the location of the springs between Vallo di Nera and Umbriano. The depth of the incision of the NE and NW-trending valleys is schematically shown in Figure 2 (left upper panel). Some of these valleys are very deep and profoundly incised. Nevertheless, they do not intercept the water table.

2.2. Recharge Estimation

Daily precipitation and mean temperature (period: 1950–2013) were collected over the study area at the station locations of 16 and 7, respectively. Data were interpolated over a 1 km2 regular grid through an ordinary kriging. Temperature data were previously detrended from the local lapse rate. Observed semivariograms were fitted through a spherical model at a monthly time step. The recharge to the aquifer was estimated, at a daily time step, over a 1 km2 square grid, using the Thornthwaite–Mather model [46,47]. Recharge was assumed to be a fraction of water surplus when soil moisture exceeds the field capacity. Soil moisture was estimated as the difference between precipitation and actual evapotranspiration, the latter computed as a fraction of the potential evapotranspiration when the soil was partially saturated. Field capacity was set to 100 mm [36]. The potential infiltration coefficients have been ascribed to each hydrostratigraphic unit on the grounds of existing technical reports on the study area [48] and on an expert judgement basis, assuming that limestone formations have higher potential infiltration coefficients than marl formations and range from 5 to 90% of the water surplus (Figure 1B). Both field capacity and potential infiltration coefficients have been calibrated to match the observed aquifer contribution to the river and springs (3.4–3.6 m3/s). The difference between water surplus and infiltration was ascribed to runoff. The validation was based on the comparison of the calculated annual infiltration with the observed aquifer discharge. The average estimated infiltration rate after calibration is 480 mm/y, corresponding to an average discharge of 3.5 m3/s. Nonstationarity in the estimated recharge was assessed through a singular spectrum analysis (SSA) for the past 60 years (1950–2012) indicating an approximately linear reduction in the recharge without evident break points. The negative trend of the recharge is 1.8 mm/y, which is significant at 90% (Figure 4). In Figure 4, the aquifer discharge from 1997 to 2012, as estimated in Section 2.1, is also shown.

3. Numerical Model Description

3.1. Layers Discretization, Boundary Conditions and Codes

The steady-state models were implemented by gradually increasing the number of layers, from a 1-layer (2D) to a 3-layer (quasi-3D) and then a 5-layer model (fully 3D) with a uniform grid spacing of 100 × 100 m (Figure 5). The top and bottom of the layers were built using the available geological data (boreholes stratigraphy, geological maps, cross sections), by means of an ordinary kriging algorithm, with a spherical semivariogram. No-flow boundary conditions were assigned to the external boundary of the models. As the upper layers do not cover the whole extent in the quasi-3D and fully 3D models, no-flow boundary conditions (inactive cells) were assigned in the areas where the anticline is eroded (eastern portion), leading to the lowermost layers to outcrop. This allows the partitioning of the recharge into different layers according to the geological setting, assigning the calculated recharge rates to each outcropping cell. Drain and river packages (DRN and RIV) were used to simulate head-dependent flux boundary conditions for the springs, the Nera River and the Precetto stream, respectively. Initial hydraulic conductivity and transmissivity values were set based on the available on-site tests and on the results of previous numerical simulations [36,48]. The 2D model features one layer only. Simulating a unique aquifer, it neglects the behavior of the two aquitards. The groundwater flow is only in the x–y direction. The hydraulic conductivity values (Kh) were set for each cell by considering the average aquifer transmissivity divided by the cell thickness. The quasi-3D model consists of three layers, representing the sub-aquifers. Aquitards are not explicitly represented; the groundwater flow in each layer is in the x–y direction and exchanges between layers are regulated through vertical conductance known as the term vertical leakage. Vertical leakage is computed by the preprocessor Groundwater Vistas (Environmental Simulations International®) based on the saturated thickness and vertical K assigned to the implicit aquitard layer [49]. Finally, the fully 3D approach allows for vertical flow in aquifers and aquitards through the explicit representation of the aquitards. The fully 3D model was set up with five layers representing the three sub-aquifers and the two aquitards. Initially, the hydraulic conductivity of the explicit aquitards was assumed to be 1/100 of the hydraulic conductivity values set to each overlaying sub-aquifer. The mean estimated recharge (3.5 m3/s, see Section 2.2) was uniformly assigned as the input recharge to the active cells of the models. Two-dimensional and quasi-3D models were run using MODFLOW 88/98 with a preconjugated gradient 2 solver (PCG2). In order to overcome dry-cell problems and reduce the model error in the fully 3D model, the Newton–Raphson formulation for MODFLOW-2005 (MODFLOW2005-NWT, [26]) was invoked. MODFLOW2005-NWT uses an alternative formulation of the GW-flow equation: the upstream weighting package (UPW) treats nonlinearities of cell drying and rewetting by using a continuous function of hydraulic head, instead of the discrete approach applied by the block-centered flow and layer property flow packages in the previous MODFLOW versions. Application of MODFLOW-NWT overcomes numerical problems by smoothing the transition from wet to dry cells and keeps all cells active [50]. MODFLOW-NWT keeps all cells active that are active at the start of the simulation. It assigns a head to unconfined cells even when the head falls below the cell bottom, allowing vertical flow in the form of recharge. However, dry cells no longer participate in horizontal aquifer flow. Use of the MODFLOW2005-NWT avoids solver instability in the presence of dry cells and diminishes the sensitivity of the solution to initial conditions.
In order to allow for a robust comparison of the 2D and quasi-3D simulations to the fully 3D model, an equivalent transmissivity was calculated for each cell for the quasi-3D and 2D models, taking account of the reduced number of layers:
T 1 + T 2 + T 3 + T 4 + T 5 = T 1 ′ + T 2 ′ + T 3 ′ = T 2 D = K 1 e 1 + K 2 e 2 + K 3 e 3 + K 4 e 4 + K 5 e 5 = K 1 ′ e 1 ′ + K 2 ′ e 2 ′ + K 3 ′ e 3 ′ = K 2 D e 2 D
  • T1, T2, T3, T4, T5 = transmissivity of the fully 3D model (subscript indicates the layer)
  • T1′, T2′, T3′ = transmissivity of the quasi-3D model (subscript indicates the layer)
  • T2D = transmissivity of the 2D model
  • Kn = permeability (subscript indicates the layer)
  • en = thickness (subscript indicates the layer)

3.2. Model Calibration

Models were calibrated in steady state through a manual trial-and-error procedure, using the available data, which comprise:
  • the discharge increment in Nera River reaches R1 to R4;
  • the discharge of Precetto stream;
  • the discharge of Lupa, Scheggino and Pacce springs (Table 1).
  • the head measured in two wells (Renari di Capriglia and Scheggino, Figure 1)
Flux and head targets used for model calibration are listed in Table 2. The topography elevation, especially in deep gorges, where the geological formations of the lower sub-aquifers are outcropping provided additional head constraints. The calibrated parameters were the horizontal and vertical hydraulic conductivities of each layer and river-bed hydraulic conductivity and thickness (Table 3 and Table 4). The streambed conductance, which enters in the computation of the groundwater–surface water interaction, is computed by the software as the ratio of riverbed hydraulic conductivity and thickness.
The calibration phase involved an iterative refining of both aquifer K values and distribution and also riverbed Kv and thickness values, on the basis of a comparison between simulation results and observations. The calibration of hydraulic conductivity in each layer was performed by gradually modifying the K values by “zones”, with each zone representing a homogeneous area of the layer. At first, the values from [36] were assigned to the sub-aquifers of the fully 3D model. During the calibration phase, the number of hydraulic conductivity zones was increased in order to refine the K distribution until the simulated values were close to the observations. After each K value variation, the model was rerun to compare results to the previous setting. Then, the pattern of K values was progressively refined, and many different K zones were added. The most difficult task was to reproduce the discharge of the river and maintain a sufficient discharge in the upper reaches. In order to reproduce this discharge distribution, a high conductivity strip was added, oriented N–S in layer 5 of the fully 3D model, which simulates faults bounding Calcare Massiccio Fm, acting as a preferential flow zone. The transmissivity distribution resulting in the fully 3D model after this calibration was transposed to the quasi-3D and 2D models following the Equation (1). Two-dimensional and quasi-3D models were further calibrated after this step in order to seek the best match using observations. In the final setting, the numbers of the calibrated conductivity zones was 3, 6 and 12 for the 2D, quasi-3D and fully 3D models, respectively (Table 3 and Figure 6).
The Nera River can be conceptualized as a bedrock stream with a limited thickness of loose alluvial sediments, up to a few meters of gravels and sands with very high permeability. In the calculation of the groundwater–surface water interaction, what really matters is the ratio between the riverbed Kv/thickness. During calibration, in order to obtain the observed high outflow to the river, this ratio was increased from 5.7 × 10−5 m/s to 7.7 × 10−4 m/s to obtain a minimal resistance to the groundwater flow from the aquifer to the river, assuming 3 m of riverbed thickness and 2.3 × 10−3 m/s for Kv (Table 4), which is consistent with the high permeability of the gravels and the presence of a few meters of riverbed above the bedrock.
Time of calculation ranges from a few seconds for the 2D model to several minutes for the fully 3D.

4. Results and Discussion

All the calibrated models give a good match between simulated and observed heads in the two only available wells (Figure 7 and Table 5). However, the 2D model was never able to match the observed distribution of river baseflow along each reach of the Nera River. In fact, assuming a single aquifer, this model is not able to feed the upstream cells of the river, and the groundwater converges toward the most downstream RIV cells (Figure 8A), which does not correspond to the observed groundwater–surface water interaction. The upstream reach R1 (Figure 1) loses water to the aquifer instead of gaining, while the downstream reach R4 gains the double of the observed discharge (Table 5).
In order to reproduce the discharge correctly partitioned along the reaches, it was necessary to construct multilayer models with consideration of system aquitards.
An important advantage of quasi-3D and fully 3D models with respect to the 2D model is that the slope and dipping of single layers can be correctly reproduced, allowing flow to be reliably simulated along the layers. In addition, the quasi-3D and fully 3D models are able to represent localized areas of preferential upward flow with a nonuniform vertical leakage array or by adding high vertical permeability zones to simulate discontinuities in the aquitards, respectively.
Further, the quasi-3D and fully 3D models allow partitioning of the recharge into different numerical layers according to the geological setting. For example, where the anticline is eroded, the recharge is assigned to the outcropping layers; thus, each of the sub-aquifers directly receive a share of the total recharge (Figure 5). In the quasi-3D and fully 3D models, the simulations show that a large part of the anticline hosts dry cells due to the high elevation of the aquifer bottom with respect to the calculated heads. When using the PCG2 solver (quasi-3D model), the dry cells are excluded from the head calculations and recharge is passed to an underlying active cell. Conversely, in the NWT upstream weighting package (fully 3D model), dry cells remain active with a calculated head value that falls below the bottom of the cell; that head is used to set up the gradient to pass recharge to an underlying cell hosting the water table or toward adjacent cells in the case of the bottom layer. In particular, the quasi-3D model performs better than the 2D, but the simulated discharge of the reach R4 is still not acceptable (Table 5). The fully 3D model shows the best performance in terms of bias (Figure 7). The hydraulic heads calculated for the three aquifers in the fully 3D model (Figure 8B) show the important role of the basal sub-aquifer (layer five) to direct the high infiltration from the mountainous recharge areas to the northern, downstream, sectors of the structure. Due to the complex geometry, sub-aquifers one and two show extended dry cells areas, and the head potentials are strongly influenced by the RIV cells.
Indeed, the use of five layers allows for a very detailed setting of the distribution of K zones which, once calibrated, should reflect a possible pattern of permeability values, making the model results consistent with the observed discharge rates and heads. Previous runs, performed with relatively high K strips in the sub-aquifers three and five, e.g., the NW–SE faults cutting the structure near Scheggino (see Figure 1 and Figure 6), were unable to drive the fluxes toward the northern part of the structure. Finally, the insertion of the high conductivity strip in the sub-aquifer three (Figure 6) ameliorated by far the calibration of the fully 3D model. Thus, the high hydraulic conductivity zones in a generalized and simplified view of the geological setting are likely to reflect fractures or faults zones. Tectonics also drive the vertical exchanges among the sub-aquifers through the aquitards. A better match was achieved considering zones of relatively high vertical permeability within the aquitards, simulating highly tectonized zones or stratigraphic gaps, for example in correspondence of reach R2 (Figure 6), which allows the groundwater to flow vertically upward from layer five to the upper layers.
However, high hydraulic conductivity zones alone were not sufficient to maintain the hydraulic gradients steep enough to feed the most upgradient reach R1. The groundwater mostly flows in the sub-aquifer three. Despite a prevalent diffuse circulation (which justifies the EPM approach), the high hydraulic conductivity strip along the anticline axis that facilitates the groundwater to flow from south to north is consistent with a discrete groundwater circulation pattern. This was confirmed by observation during the perforation of the Renari well for the Calcare Massiccio aquifer, which showed a prevalent circulation through fractures. Table 5 lists the results of the three models after calibration. Figure 7 shows the bias %.
Ultimately, the calibration process points to the presence of structural elements in the flow system that are not readily observable by other means. In this sense the stepwise process is not meant to produce a unique representation of the subsurface, but rather points to the existence of subsurface features that seem to be controlling flow to the river according to the fully 3D model, but which are otherwise very difficult to characterize by field work, i.e., the examination of outcrops or geophysical prospections. Such a model can be considered as an interpretive model [44] in the sense of a screening model that helps the modeler to develop an initial understanding of a groundwater system and/or test hypotheses about the system. Calibration through manual trial-and-error proved that a simple 2D model is not able to match field observations and thus is not an acceptable model for the site. Further, it provided insights into how parameter changes in different areas of the model could correspond to unknown subsurface features [51]. For example, the high permeability strip in the deepest layer in the fully 3D model could suggest the existence of a still unknown, developed karstic system. This manual process, although seemingly unable to find a unique solution, might be the only way forward given the impossibility of quantifying the “true uncertainty” in natural systems and the error inevitably associated with a complex structural model with few data supporting observations [52]. As stated by Fienen at al. (2009) [53], parsimonious models should avoid unnecessary or unsupported complexity while accurately delineating flow paths given our state of knowledge about the field setting.

5. Conclusions

Modelling a complexly folded and faulted hydrogeological system requires a stepwise procedure to reach the best solution with a parsimonious model setting. The inherent complexity of the conceptual model was reproduced in the numerical model by progressively adding elements to the model grid such as new layers from the 2D to the fully 3D model and hydraulic conductivity zones, and by calibrating the steady-state solution by manual trial and error. Furthermore, gradually increasing the model complexity can be a useful approach for simulating groundwater flow when few head data are available as it is common in karst aquifers. Due to the prevalent diffuse circulation, an EPM approach was used; however, the calibrated setting of hydraulic conductivity zones suggests a discrete groundwater circulation pattern, which was successfully simulated by adding a high permeability longitudinal strip. The calibrated pattern of K zones both for sub-aquifers and aquitards is likely to reflect the structural and stratigraphic setting. The quasi-3D and the fully 3D models both allow for recharge partitioning into the sub-aquifers. The vertical exchanges among the sub-aquifers are regulated by leakage coefficient or aquitards parametrization, respectively. The higher number of layers compared to the 2D model allows the 3D simulations to drive the groundwater flow towards the different parts of the river. However, the fully 3D model best matches the observed flow distribution at the different reaches along the river, simulating reliable flow paths and recharge partitioning into layers. The Newton–Raphson formulation of MODFLOW2005 is required to achieve convergence and reduce model error mainly due to cells drying and rewetting and producing numerical instabilities. The numerical model demonstrated the major impact of folded and faulted geological structures on controlling the flow dynamics in terms of flow direction, water heads and spatial distribution of the outflows to the river and springs. The stepwise process of model construction and calibration, even with a limited number of head and flux targets, points to the presence of structural elements in the subsurface that otherwise can escape observation in field studies of the terrain.

Author Contributions

Conceptualization, E.P. and C.D.S.; methodology, E.P., C.D.S., N.G. and E.R.; software, E.P., C.D.S., A.B.P. and N.G.; validation, E.P., C.D.S. and E.R.; formal analysis, E.P., C.D.S. and E.R.; investigation, E.P. and C.D.S.; resources, E.P.; data curation, E.P. and C.D.S.; writing—original draft preparation, E.P. and C.D.S.; writing—review and editing, E.P., C.D.S. and E.R.; visualization, C.D.S., A.B.P., E.R. and N.G.; supervision, E.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data of springs discharge have been provided by the Regional Agency for Environmental Protection (ARPA Umbria). Precipitation, temperature and river discharge data have been provided by the Hydrographic Service of Umbria (Ufficio Idrografico della Regione Umbria).

Acknowledgments

The authors are deeply grateful to Daniel Feinstein and Rundy Hunt (USGS) for their useful suggestions and review of the manuscript. We also wish to thank Reviewer 1 for his/her interesting comments and constructive criticisms.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Ford, D.C.; Williams, P.W. Karst Hydrogeology and Geomorphology; John Wiley & Sons: Chichester, UK; Hoboken, NJ, USA, 2007. [Google Scholar]
  2. Worthington, S.; Ford, D. Self-Organized Permeability in Carbonate Aquifers. Ground Water 2009, 47, 326–336. [Google Scholar] [CrossRef] [PubMed]
  3. Worthington, S.R.H. Management of Carbonate Aquifers. In Karst Management; Van Beynen, P.E., Ed.; Department Geography, University of South Florida: Tampa, FL, USA, 2011; Volume XII, pp. 243–261. [Google Scholar] [CrossRef]
  4. Dragoni, W.; Mottola, A.; Cambi, C. Modeling the effects of pumping wells in spring management: The case of Scirca spring (central Apennines, Italy). J. Hydrol. 2013, 493, 115–123. [Google Scholar] [CrossRef]
  5. Fiorillo, F.; Guadagno, F.M. Karst Spring Discharges Analysis in Relation to Drought Periods, Using the SPI. Water Resour. Manag. 2010, 24, 1867–1884. [Google Scholar] [CrossRef]
  6. Romano, E.; Preziosi, E. Precipitation pattern analysis in the Tiber River basin (central Italy) using standardized indices. Int. J. Climatol. 2013, 33, 1781–1792. [Google Scholar] [CrossRef]
  7. Citrini, A.; Camera, C.; Beretta, G.P. Nossana Spring (Northern Italy) under Climate Change: Projections of Future Discharge Rates and Water Availability. Water 2020, 12, 387. [Google Scholar] [CrossRef]
  8. Leone, G.; Pagnozzi, M.; Catani, V.; Ventafridda, G.; Esposito, L.; Fiorillo, F. A hundred years of Caposele spring discharge measurements: Trends and statistics for understanding water resource availability under climate change. Stoch. Hydrol. Hydraul. 2020, 35, 345–370. [Google Scholar] [CrossRef]
  9. Romano, E.; Petrangeli, A.B.; Salerno, F.; Guyennon, N. Do recent meteorological drought events in central Italy result from long-term trend or increasing variability? Int. J. Clim. 2021, 42, 4111–4128. [Google Scholar] [CrossRef]
  10. Romano, E.; Petrangeli, A.B.; Preziosi, E. Spati al and Time Analysis of Rainfall in the Tiber River Basin (Central Italy) in relation to Discharge Measurements (1920-2010). Procedia Environ. Sci. 2011, 7, 258–263. [Google Scholar] [CrossRef]
  11. IPCC 2022—Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change; Pörtner, H.-O.; Roberts, D.C.; Tignor, M.; Poloczanska, E.S.; Mintenbeck, K.; Alegría, A.; Craig, M.; Langsdorf, S.; Löschke, S.; Möller, V.; et al. (Eds.) Cambridge University Press: Cambridge, UK, 2022. [Google Scholar]
  12. Preziosi, E.; Romano, E. Are large karstic springs good indicators for Climate Change effects on groundwater? In Proceedings of the European Geosciences Union General Assembly 2013, Vienna, Austria, 7–12 April 2013. Geophys. Res. Abstr. 2013, 15, 9367. [Google Scholar]
  13. Romano, E.; Del Bon, A.; Petrangeli, A.; Preziosi, E. Generating synthetic time series of springs discharge in relation to standardized precipitation indices. Case study in Central Italy. J. Hydrol. 2013, 507, 86–99. [Google Scholar] [CrossRef]
  14. Band, S.S.; Heggy, E.; Bateni, S.M.; Karami, H.; Rabiee, M.; Samadianfard, S.; Chau, K.-W.; Mosavi, A. Groundwater level prediction in arid areas using wavelet analysis and Gaussian process regression. Eng. Appl. Comput. Fluid Mech. 2021, 15, 1147–1158. [Google Scholar] [CrossRef]
  15. Mallick, J.; Almesfer, M.K.; Alsubih, M.; Talukdar, S.; Ahmed, M.; Kahla, N.B. Developing a new method for future groundwater potentiality mapping under climate change in Bisha watershed, Saudi Arabia. Geocarto Int. 2022, 1, 1–33. [Google Scholar] [CrossRef]
  16. Snoussi, M.; Jerbi, H.; Tarhouni, J. Integrated Groundwater Flow Modeling for Managing a Complex Alluvial Aquifer Case of Study Mio-Plio-Quaternary Plain of Kairouan (Central Tunisia). Water 2022, 14, 668. [Google Scholar] [CrossRef]
  17. European Community. Directive 2000/60/EC of the European Parliament and of the Council of 23 October 2000 Establishing a Framework for Community Action in the Field of Water Policy. 2000. Available online: https://eur-lex.europa.eu/LexUriServ/LexUriServ.do?uri=CELEX:32000L0060:EN:HTML (accessed on 5 July 2022).
  18. Gattinoni, P.; Francani, V. Depletion risk assessment of the Nossana Spring (Bergamo, Italy) based on the stochastic modeling of recharge. Appl. Hydrogeol. 2010, 18, 325–337. [Google Scholar] [CrossRef]
  19. Levison, J.; Larocque, M.; Ouellet, M. Modeling low-flow bedrock springs providing ecological habitats with climate change scenarios. J. Hydrol. 2014, 515, 16–28. [Google Scholar] [CrossRef]
  20. Reimann, T.; Giese, M.; Geyer, T.; Liedl, R.; Maréchal, J.C.; Shoemaker, W.B. Representation of water abstraction from a karst conduit with numerical discrete-continuum models. Hydrol. Earth Syst. Sci. 2014, 18, 227–241. [Google Scholar] [CrossRef]
  21. Scanlon, B.R.; Mace, R.E.; Barrett, M.E.; Smith, B. Can we simulate regional groundwater flow in a karst system using equivalent porous media models? Case study, Barton Springs Edwards aquifer, USA. J. Hydrol. 2003, 276, 137–158. [Google Scholar] [CrossRef]
  22. Abusaada, M.; Sauter, M. Studying the Flow Dynamics of a Karst Aquifer System with an Equivalent Porous Medium Model. Ground Water 2012, 51, 641–650. [Google Scholar] [CrossRef]
  23. Ben-Itzhak, L.; Gvirtzman, H. Groundwater flow along and across structural folding: An example from the Judean Desert, Israel. J. Hydrol. 2005, 312, 51–69. [Google Scholar] [CrossRef]
  24. Dafny, E.; Burg, A.; Gvirtzman, H. Effects of Karst and geological structure on groundwater flow: The case of Yarqon-Taninim Aquifer, Israel. J. Hydrol. 2010, 389, 260–275. [Google Scholar] [CrossRef]
  25. Valota, G.; Giudici, M.; Parravicini, G.; Ponzini, G.; Romano, E. Is the forward problem of ground water hydrology always well posed? Ground Water 2002, 40, 500–508. [Google Scholar] [CrossRef] [PubMed]
  26. Niswonger, R.G.; Panday, S.; Ibaraki, M. MODFLOW-NWT, A Newton formulation for MODFLOW-2005. US Geol. Surv. Tech. Methods 2011, 6, 44. [Google Scholar]
  27. Mantoglou, A. On optimal model complexity in inverse modeling of heterogeneous aquifers. J. Hydraul. Res. 2005, 43, 574–583. [Google Scholar] [CrossRef]
  28. Hill, M.C. The practical use of simplicity in developing groundwater models. Ground Water 2006, 44, 775–781. [Google Scholar] [CrossRef] [PubMed]
  29. Moore, C.; Wöhling, T.; Doherty, J. Efficient regularization and uncertainty analysis using a global optimization methodology. Water Resour. Res. 2010, 46, W08527. [Google Scholar] [CrossRef]
  30. Hill, M.C. Methods and Guidelines for Effective Model Calibration. In U.S. Geological Survey, Water-Resources Investigations Report 98–4005; U.S. Geological Survey: Denver, CO, USA, 1998; p. 98. Available online: https://water.usgs.gov/nrp/gwsoftware/modflow2000/WRIR98-4005.pdf (accessed on 10 August 2022).
  31. Barchi, M.R.; Lemmi, M. Geology of the area of Monte Coscerno-Monte di Civitella (South-eastern Umbria). Boll. Soc. Geol. Ital. 1996, 115, 601–624. [Google Scholar]
  32. Tavarnelli, E. Analisi geometrica e cinematica dei sovrascorrimenti compresi fra la Valnerina e la conca di Rieti (Appennino Umbro-Marchigiano-Sabino). Boll. Soc. Geol. Ital. 1996, 113, 249–259. [Google Scholar]
  33. Pace, P.; Calamita, F.; Tavarnelli, E. Brittle-Ductile Shear Zones along Inversion-Related Frontal and Oblique Thrust Ramps: Insights from the Central-Northern Apennines Curved Thrust System (Italy). In Ductile Shear Zones: From Micro- to Macro-Scales; Mukherjee, S., Mulchrone, K.F., Eds.; John Wiley & Sons: Chichester, UK, 2015; Chapter 8; pp. 111–127. [Google Scholar]
  34. Calamita, F.; Deiana, G. The arcuate shape of the Umbria-Marche-Sabina Apennines (central Italy). Tectonophysics 1988, 146, 139–147. [Google Scholar] [CrossRef]
  35. Boni, C.; Preziosi, E. Una possibile simulazione numerica dell’acquifero basale di M.Coscerno—M.Aspra (bacino del fiume Nera). Geol. Appl. E Idrogeol. 1993, 28, 131–140. [Google Scholar]
  36. Preziosi, E. Simulazioni Numeriche di Acquiferi Carbonatici in Aree Corrugate: Applicazioni al Sistema Idrogeologico Della Valnerina (Italia Centrale). In Quad IRSA-CNR; Istituto di Ricerca Sulle Acque-CNR: Rome, Italy, 2007; Volume 125, p. 225. ISSN 0390-6329. [Google Scholar]
  37. Mastrorillo, L.; Baldoni, T.; Banzato, F.; Boscherini, A.; Cascone, D.; Checcucci, R.; Petitta, M.; Boni, C. Quantitative hydrogeological analysis of the carbonate domain of the Umbria Region (Central Italy). Ital. J. Eng. Geol. Environ. 2009, 1, 137–155. [Google Scholar] [CrossRef]
  38. Mastrorillo, L.; Petitta, M. Hydrogeological conceptual model of the upper Chienti River Basin aquifers (Umbria-Marche Apennines). Ital. J. Geosci. 2014, 133, 396–408. [Google Scholar] [CrossRef]
  39. Di Matteo, L.; Dragoni, W.; Ulderico, V.; Latini, M.; Spinsanti, R. Risorse idriche sotterranee e loro gestione: Il caso dell’ATO2 Umbria (Umbria meridionale). Acque Sotter. 2005, 96, 9–21. [Google Scholar]
  40. Damiani, A.V.; Baldanza, A.; Barchi, M.R.; Boscherini, A.; Checcucci, R.; Decandia, F.A.; Felicioni, G.; Lemmi, M.; Luchetti, L.; Motti, A.; et al. Geological Map of Italy. Scale 1:50,000. Scheet 563 “Spoleto”; ISPRA-Servizio Geologico Nazionale: Roma, Italy, 2008. [Google Scholar]
  41. ARPA Umbria. Available online: https://apps.arpa.umbria.it/acqua/contenuto/monitoraggio-in-continuo-acque-sotterranee (accessed on 13 June 2022).
  42. Preziosi, E.; Romano, E. From hydrostructural analysis to mathematical modelling of regional aquifers. Ital. J. Eng. Geol. Environ. 2009, 1, 183–198. [Google Scholar] [CrossRef]
  43. ARPA Umbria. Il Sistema Informativo di Arpa Umbria e la Conoscenza Ambientale, a cura di Mauro Emiliano, Quaderni Arpa Umbria. 2004. Available online: https://www.arpa.umbria.it/pagine/il-sistema-informativo-di-arpa-umbria-e-la-conosce (accessed on 5 July 2022).
  44. White, W.B. Geomorphology and Hydrology of Karst Terrains; Oxford Universuty Press: New York, NY, USA, 1988; p. 464. [Google Scholar]
  45. Vigna, B.; Banzato, C. The hydrogeology of high-mountain carbonate areas: An example of some Alpine systems in southern Piedmont (Italy). Environ. Earth Sci. 2015, 74, 267–280. [Google Scholar] [CrossRef]
  46. Thornthwaite, C.W.; Mather, J.R. The water balance. Climatology 1955, 8, 5–86. [Google Scholar]
  47. Thornthwaite, C.W.; Mather, J.R. Instructions and Tables for Computing Potential Evapotranspiration and the Water Balance. Centerton 1957, 10, 185–310. [Google Scholar]
  48. Umbriadue Servizi Idrici S.C.A R.L.; Technologies for Water Services S.p.A. Acquedotto Scheggino-Pentima. In Progetto Definitivo, Technical Report; Servizio Idrico Integrato: Terni, Italy, 2006. [Google Scholar]
  49. Rambaugh, J.O.; Rambaugh, D.B. Guide to Using Groundwater Vistas, Version 7, Environmental Simulation. 2017. Available online: https://d3pcsg2wjq9izr.cloudfront.net/files/4134/download/730323/4-gv7manual.pdf (accessed on 5 July 2022).
  50. Feinstein, D.T.; Fienen, M.N.; Kennedy, J.L. Development and Application of a Groundwater/Surface-Water Flow Model Using MODFLOW-NWT for the Upper Fox River Basin, Southeastern Wisconsin. In U.S. Geological Survey Scientific Investigations Report; CreateSpace Independent Publishing Platform: Scotts Valley, CA, USA, 2012. [Google Scholar]
  51. Anderson, M.P.; Woessner, W.W.; Hunt, R.J. Applied Groundwater Modeling: Simulation of Flow and Advective Transport, 2nd ed.; Academic Press, Inc.: Cambridge, MA, USA, 2015; p. 564. ISBN 9780120581030. [Google Scholar]
  52. Hunt, R.J. Applied Uncertainty. Groundwater 2017, 55, 771–772. [Google Scholar] [CrossRef]
  53. Fienen, M.F.; Hunt, R.J.; Krabbenhoft, D.P.; Clemo, T. Obtaining parsimonious hydraulic conductivity fields using head and transport observations: A Bayesian Geostatistical parameter estimation approach. Water Resour. Res. 2009, 45, W08405. [Google Scholar] [CrossRef]
Water 14 02475 g001 550
Figure 1. (A): Hydrogeological setting, conceptual model. (B): distribution map of the potential infiltration coefficients. Tectonic elements (faults, thrust, anticlinal axis) from [40].
Figure 1. (A): Hydrogeological setting, conceptual model. (B): distribution map of the potential infiltration coefficients. Tectonic elements (faults, thrust, anticlinal axis) from [40].
Water 14 02475 g001
Water 14 02475 g002 550
Figure 2. Topographic map of the study area. Elevation in meters above sea level.
Figure 2. Topographic map of the study area. Elevation in meters above sea level.
Water 14 02475 g002
Water 14 02475 g003 550
Figure 3. Schematic cross section and translation into model grid for the fully 3D model. Trace of cross section and legend are in Figure 1. (Top panel): schematic hydrogeologic cross section. The blue filling at the top and central panels represents the saturated zone; numbers 1, 2 and 3 refer to the sub-aquifers. Dotted lines in central panel represent the potentiometric levels of confined sub-aquifers. (Bottom panel): model grid for the fully 3D model; light blue arrows: recharge; blue arrow: discharge to the RIV cells.
Figure 3. Schematic cross section and translation into model grid for the fully 3D model. Trace of cross section and legend are in Figure 1. (Top panel): schematic hydrogeologic cross section. The blue filling at the top and central panels represents the saturated zone; numbers 1, 2 and 3 refer to the sub-aquifers. Dotted lines in central panel represent the potentiometric levels of confined sub-aquifers. (Bottom panel): model grid for the fully 3D model; light blue arrows: recharge; blue arrow: discharge to the RIV cells.
Water 14 02475 g003
Water 14 02475 g004 550
Figure 4. Estimated recharge to the aquifer (purple graph), with respect to the cumulated precipitation (histogram) and the observed aquifer discharge (blue graph). Trend (black dashed line) and Sen’s slope refer to the calculated recharge.
Figure 4. Estimated recharge to the aquifer (purple graph), with respect to the cumulated precipitation (histogram) and the observed aquifer discharge (blue graph). Trend (black dashed line) and Sen’s slope refer to the calculated recharge.
Water 14 02475 g004
Water 14 02475 g005 550
Figure 5. Sketch of the three steady-state models implemented with 2D, quasi-3D and fully 3D settings. Only active cells are shown.
Figure 5. Sketch of the three steady-state models implemented with 2D, quasi-3D and fully 3D settings. Only active cells are shown.
Water 14 02475 g005
Water 14 02475 g006 550
Figure 6. Distribution of K zones for fully 3D calibrated model. R1, R2, R3 and R4 indicates the river reaches (see also Figure 1).
Figure 6. Distribution of K zones for fully 3D calibrated model. R1, R2, R3 and R4 indicates the river reaches (see also Figure 1).
Water 14 02475 g006
Water 14 02475 g007 550
Figure 7. Target bias on the discharge (%).
Figure 7. Target bias on the discharge (%).
Water 14 02475 g007
Water 14 02475 g008 550
Figure 8. Simulated head for the 2D model (A) and for the layers 1, 3 and 5 in the fully 3D model (B). Model cross sections at the top of figures (row 162, column 23). Approximate locations of Scheggino well (S) and Renari well (R) are shown.
Figure 8. Simulated head for the 2D model (A) and for the layers 1, 3 and 5 in the fully 3D model (B). Model cross sections at the top of figures (row 162, column 23). Approximate locations of Scheggino well (S) and Renari well (R) are shown.
Water 14 02475 g008
Table
Table 1. Gauging stations and springs in the study area.
Table 1. Gauging stations and springs in the study area.
Site NameAltitude (m a.s.l.)Average Discharge (L/s)Reference PeriodReference
Borgo Cerreto34524511991–1993[36]
Vallo di Nera29328871991–1993[36]
Scheggino27540071991–1993[36]
Ceselli26551171991–1993[36]
Umbriano24257171991–1993[36]
T.Orsina21160171991–1993[36]
Scheggino spring2761901991–1993(accounted for in the Ceselli gauging site) [36]
Lupa spring3661201997–2012[41]
Pacce spring480602000–2001Well field excluded. [43]
Precetto stream3251201991–1993[36]
Borgo Cerreto34524511991–1993[36]
Vallo di Nera29328871991–1993[36]
Table
Table 2. Flux and head targets used for model calibration.
Table 2. Flux and head targets used for model calibration.
LayerHead TargetFlux Target
1 (sub-aquifer 1) Scaglia unitn.a.Nera River reach R1, Lupa spring, Scheggino spring, Precetto stream
3 (sub-aquifer 2), Maiolica unitScheggino wellNera River reach R2, R3 and R4, Pacce spring
5 (sub-aquifer 3), Calcare massiccio and CorniolaRenari di Capriglia welln.a.
Table
Table 3. Range of calibrated Kh and Kz values.
Table 3. Range of calibrated Kh and Kz values.
Model SettingN of LayersN of Active CellsN of Calibrated K ZonesRange of Kh Values, m/sRange of Kz Values, m/s
2D123,24831.1 × 10−5/2.8 × 10−51.15 × 10−6/2.8 × 10−6
Quasi-3D350,95960.1 × 10−6/4.6 × 10−51.1 × 10−7/2.8 × 10−6
Fully 3D577,212121.1 × 10−6/8.1 × 10−5 (high K strip: 5.7 × 10−4 m/s)1.1 × 10−7/8.1 × 10−5 (high K strip: 6.3 × 10−4 m/s)
Table
Table 4. Initial and calibrated values of riverbed hydraulic conductivity and thickness.
Table 4. Initial and calibrated values of riverbed hydraulic conductivity and thickness.
River ReachRiverbed Kv (m/s)Riverbed Thickness (m)
InitialCalibratedInitialCalibrated
1–4 (Nera River)5.7 × 10−52.3 × 10−313.0
5 (Precetto stream)5.7 × 10−52.3 × 10−313.0
Table
Table 5. Results of the three models and comparison with the observed values. The bias % is shown.
Table 5. Results of the three models and comparison with the observed values. The bias % is shown.
ElementObserved2DQuasi-3DFully 3D
Well Scheggino, m a.s.l.273271.51275.57273.26
Well Renari, m a.s.l.301269.88312.52305.00
Nera R., R1, l/s436−322.00351.92238.56
Nera R., R2, l/s11201557.82812.62778.68
Nera R., R3, l/s1110507.18751.171445.15
Nera R., R4, l/s6091173.951760.49669.35
springs, l/s390587.99195.44286.21
Precetto s., l/s12060.40117.71103.29
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Preziosi, E.; Guyennon, N.; Petrangeli, A.B.; Romano, E.; Di Salvo, C. A Stepwise Modelling Approach to Identifying Structural Features That Control Groundwater Flow in a Folded Carbonate Aquifer System. Water 2022, 14, 2475. https://doi.org/10.3390/w14162475

AMA Style

Preziosi E, Guyennon N, Petrangeli AB, Romano E, Di Salvo C. A Stepwise Modelling Approach to Identifying Structural Features That Control Groundwater Flow in a Folded Carbonate Aquifer System. Water. 2022; 14(16):2475. https://doi.org/10.3390/w14162475

Chicago/Turabian Style

Preziosi, Elisabetta, Nicolas Guyennon, Anna Bruna Petrangeli, Emanuele Romano, and Cristina Di Salvo. 2022. "A Stepwise Modelling Approach to Identifying Structural Features That Control Groundwater Flow in a Folded Carbonate Aquifer System" Water 14, no. 16: 2475. https://doi.org/10.3390/w14162475

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop