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

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 km²) 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 m³/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. Gauging stations and springs in the study area.

Site Name	Altitude (m a.s.l.)	Average Discharge (L/s)	Reference Period	Reference
Borgo Cerreto	345	2451	1991–1993	[36]
Vallo di Nera	293	2887	1991–1993	[36]
Scheggino	275	4007	1991–1993	[36]
Ceselli	265	5117	1991–1993	[36]
Umbriano	242	5717	1991–1993	[36]
T.Orsina	211	6017	1991–1993	[36]
Scheggino spring	276	190	1991–1993	(accounted for in the Ceselli gauging site) [36]
Lupa spring	366	120	1997–2012	[41]
Pacce spring	480	60	2000–2001	Well field excluded. [43]
Precetto stream	325	120	1991–1993	[36]
Borgo Cerreto	345	2451	1991–1993	[36]
Vallo di Nera	293	2887	1991–1993	[36]

).

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 m³/s, ranging from 0.9 to 7.46 m³/s (Figure 4). The total aquifer discharge was estimated at about 3.4–3.6 m³/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 km² 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 km² 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 m³/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 m³/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 (K_h) 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 m³/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:

$$\begin{gathered} T_1+T_2+T_3+T_4+T_5=T_{1^{\prime }}+T_{2^{\prime }}+T_{3^{\prime }}=T_{2D}=K_1{e}_1+K_2{e}_2+K_3{e}_3+K_4{e}_4+K_5{e}_5 \\ =K_{1^{\prime }}{e}_{1^{\prime }}+K_{2^{\prime }}{e}_{2^{\prime }}+K_{3^{\prime }}{e}_{3^{\prime }}=K_{2D}{e}_{2D} \end{gathered}$$

(1)

• T₁, T₂, T₃, T₄, T₅ = transmissivity of the fully 3D model (subscript indicates the layer)
• T_1′, T_2′, T_3′ = transmissivity of the quasi-3D model (subscript indicates the layer)
• T_2D = transmissivity of the 2D model
• K_n = permeability (subscript indicates the layer)
• e_n = 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. Flux and head targets used for model calibration.

Layer	Head Target	Flux Target
1 (sub-aquifer 1) Scaglia unit	n.a.	Nera River reach R1, Lupa spring, Scheggino spring, Precetto stream
3 (sub-aquifer 2), Maiolica unit	Scheggino well	Nera River reach R2, R3 and R4, Pacce spring
5 (sub-aquifer 3), Calcare massiccio and Corniola	Renari di Capriglia well	n.a.

. 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. Range of calibrated Kh and Kz values.

Model Setting	N of Layers	N of Active Cells	N of Calibrated K Zones	Range of Kh Values, m/s	Range of Kz Values, m/s
2D	1	23,248	3	1.1 × 10−5/2.8 × 10−5	1.15 × 10−6/2.8 × 10−6
Quasi-3D	3	50,959	6	0.1 × 10−6/4.6 × 10−5	1.1 × 10−7/2.8 × 10−6
Fully 3D	5	77,212	12	1.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)

and

Table 4. Initial and calibrated values of riverbed hydraulic conductivity and thickness.

River Reach	Riverbed Kv (m/s)		Riverbed Thickness (m)
	Initial	Calibrated	Initial	Calibrated
1–4 (Nera River)	5.7 × 10−5	2.3 × 10−3	1	3.0
5 (Precetto stream)	5.7 × 10−5	2.3 × 10−3	1	3.0

). 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⁻⁵ m/s to 7.7 × 10⁻⁴ 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⁻³ 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. Results of the three models and comparison with the observed values. The bias % is shown.

Element	Observed	2D	Quasi-3D	Fully 3D
Well Scheggino, m a.s.l.	273	271.51	275.57	273.26
Well Renari, m a.s.l.	301	269.88	312.52	305.00
Nera R., R1, l/s	436	−322.00	351.92	238.56
Nera R., R2, l/s	1120	1557.82	812.62	778.68
Nera R., R3, l/s	1110	507.18	751.17	1445.15
Nera R., R4, l/s	609	1173.95	1760.49	669.35
springs, l/s	390	587.99	195.44	286.21
Precetto s., l/s	120	60.40	117.71	103.29

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