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

: This paper concerns a stepwise modelling procedure for groundwater ﬂow 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 ﬁve layers (fully 3D), using an equivalent porous medium approach, in order to ﬁnd 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 outﬂow to the river baseﬂow. The numerical model demonstrated the major impact of folded and faulted geological structures on controlling the ﬂow dynamics in terms of ﬂow direction, water heads and the spatial distribution of the outﬂows to the river and springs.


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.

Geological and Hydrogeological Setting
The case study aquifer (Monte Coscerno, 230 km 2 ) 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 (Figures 1 and 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 3 /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 m 3 /s, ranging from 0.9 to 7.46 m 3 /s ( Figure 4). The total aquifer discharge was estimated at about 3.4-3.6 m 3 /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.     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 3 /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).   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.

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

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 2 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 2 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 3 /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 3 /s. Nonstationarity in the estimated recharge was assessed through a singular spectrum analysis (SSA) for the past 60 years  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.

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 3 /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: wet to dry cells and keeps all cells active [50]. MODFLOW-NWT keeps all cells active th are active at the start of the simulation. It assigns a head to unconfined cells even whe the head falls below the cell bottom, allowing vertical flow in the form of recharg However, dry cells no longer participate in horizontal aquifer flow. Use of th MODFLOW2005-NWT avoids solver instability in the presence of dry cells an diminishes the sensitivity of the solution to initial conditions.

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 subaquifers 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 (Tables 3 and 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. meability. In the calculation of the groundwater-surface water interaction, what matters is the ratio between the riverbed Kv/thickness. During calibration, in order tain the observed high outflow to the river, this ratio was increased from 5.7 × 10 −5 7.7 × 10 −4 m/s to obtain a minimal resistance to the groundwater flow from the aqu the river, assuming 3 m of riverbed thickness and 2.3 × 10 −3 m/s for Kv (Table 4), w consistent with the high permeability of the gravels and the presence of a few me riverbed above the bedrock.
Time of calculation ranges from a few seconds for the 2D model to several m for the fully 3D. Figure 6. Distribution of K zone s for fully 3D calibrated model. R1, R2, R3 and R4 indicate s th re aches (see also Figure 1).

Results and Discussion
All the calibrated models give a good match between simulated and observed in the two only available wells (Figure 7 and Table 5). However, the 2D model was able to match the observed distribution of river baseflow along each reach of the River. In fact, assuming a single aquifer, this model is not able to feed the upstream of the river, and the groundwater converges toward the most downstream RIV cell Time of calculation ranges from a few seconds for the 2D model to several minutes for the fully 3D.

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).
Wa ter 2022, 14, x FOR PEER REVIEW 12 of 17 ure 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 subaquifers 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 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, subaquifers 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 Figures 1 and 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.

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