Modeling a Metamorphic Aquifer through a Hydro-Geophysical Approach: The Gap between Field Data and System Complexity

: The productivity of metamorphic aquifers is generally lower than that of the more common alluvial and carbonates ones. However, in some Mediterranean areas, such as the Calabria region (Italy), water scarcity combined with the presence of extensive metamorphic water bodies requires the development of further studies to characterize the hydrodynamic properties of these groundwater systems in order to achieve their sustainable exploitation. The interest in this goal becomes even greater if climate change effects are considered. The purpose of this study was to provide the geological-structural and hydrogeological numerical modeling of a metamorphic aquifer, using direct and indirect data measurement, in a large area of the Sila Piccola in Calabria. The hydrodynamic characterization of the crystalline-metamorphic aquifer, constituted by granite and metamorphic rocks, is extremely complex. The MODFLOW-2005 groundwater model was used to simulate ﬂow phenomena in the aquifer, obtaining hydraulic conductivity values of 2.7 × 10 − 6 m/s, which turned out to be two orders of magnitude higher than those obtained from the interpretation of the slug-tests performed in the study area. The mathematical model was also able to estimate the presence of a lateral recharge from a neighboring deep aquifer providing a signiﬁcant water supply to the system under investigation.


Introduction
Groundwater in altered and fractured metamorphic aquifers represents the only water resource for many areas, especially in arid and semi-arid regions [1,2].These aquifers are generally characterized by low water productivity due to their low permeability and porosity.However, the presence of water in these formations depends on various factors such as lithology, tectonics, geodynamics, and climatic conditions.Many studies have shown that their presence is inside the surface mantle and in fissured rocks where the groundwater flow takes place.Moreover, geochemical alteration, tectonic activity, and lithostatic decompression processes are the cause of the generation of cracks inside metamorphic rocks [3][4][5][6].Geodynamic and geomorphic processes affecting the presence of groundwater also have an important influence on the hydraulic conductivity of these geological formations.In crystalline rocks, for example, the hydraulic conductivity can vary up to 12 orders of magnitude [7][8][9][10].
The estimate of these parameters, at the aquifer scale, is not the only problem to face for a correct modeling of groundwater flow and transport processes, but also the reconstruction of their spatial pattern is not an easy task.Several methods for the regionalization of hydraulic conductivity or transmissivity, at the aquifer scale, have been proposed, but most of them have been applied to sedimentary or alluvial aquifers.These methods use extensive hydrodynamic parameter datasets, deduced from pumping or slug tests, to obtain their spatial distribution through geostatistical approaches [11][12][13][14][15][16][17].Geostatistical techniques can be combined with geological facies information to reconstruct the aquifer, as a multiple continuum or a composite medium, providing promising results for the upscaling of hydrodynamic parameters [18][19][20][21][22][23][24][25].Further methods based on hydro-geophysical approaches were developed, e.g., [26][27][28][29][30], leading to transmissivity regionalization by combining them with geostatistical techniques, e.g., [31,32].Such a setup, however, requires thorough field investigations [33].
The characterization of the hydrodynamic properties in hard-rock aquifers, i.e., made up of granite and metamorphic rocks, is a complex problem due to the strong heterogeneity of the system.A network of secondary connections, of different order and degree, takes place among the main fractures, generating a continuous variation in the hydrodynamic properties at different scales and defining the modulus and direction of the groundwater flow, e.g., [34][35][36][37].The field data amount that would be necessary to properly define such systems is difficult to achieve, so the continuum approach is often used in order to deal with large-scale hydrodynamic parameter spatial distributions [38,39].Although several studies have concluded that the approximation of a fractured medium into a porous one is only valid in densely fractured aquifers [40][41][42], the scientific debate on this topic is still very open, especially regarding the scale at which this approximation is still suitable [43][44][45][46][47].
In Italy, large areas of the country are covered by metamorphic rocks, and yet the hydrogeological properties of these formations are not well known.Some of the few studies, carried out in the national territory through injection tests, clearly showed the extreme heterogeneity of the rock mass, highlighted by the wide range of values determined for the hydraulic conductivity (10 −8 ÷ 10 −6 m/s) [3,48].The knowledge of this parameter is of fundamental importance given the extreme variability and heterogeneity of the lithologies with which an aquifer may be constituted.Figure 1 shows representative values of hydraulic conductivity for selected unconsolidated sedimentary materials, and sedimentary and crystalline rocks [49,50].In some Mediterranean areas such as in Calabria region, the scarcity of water resources and the considerable extension of metamorphic aquifers (39% of the regional area) would require further studies either on their hydrodynamic properties and their hydraulic behavior in order to achieve their sustainable exploitation.
Interest in these metamorphic aquifers becomes even greater if climate change effects are considered.Many climatic scenarios relating to the Mediterranean area foresee a reduction in rainfall and an increase in storm water events; this will lead to a lower infiltration recharge for aquifers and a higher hydrogeological hazard for the territory [51].This is the context of our work focusing on the study of a metamorphic aquifer in the Calabrian Arc where a reduction in the average water availability, from one side, and a high risk of landslides in the wettest periods, on the other, represent a major challenge to deal with [52].
The purposes of the study are: (i) to contextualize the geological and structural setting and the gravitational phenomena occurring in a large area of Sila Piccola, Calabria; (ii) to recognize and characterize the geometries of the metamorphic aquifer, located in that zone, using field data derived from springs, wells, and boreholes; (iii) to define the geometricstructural relationships of the geological and lithological units.

Materials and Methods
The basin including the slope on the inhabited center of Gimigliano was the object of in-depth study through the realization of continuous-core boreholes.Inclinometers and piezometers were installed at different depths.The obtained data were continuously processed for a significant time interval.Mass movements of variable size and depth were mapped.Particular attention was paid to large landslides, together with other geomorphological processes, which have conditioned the morphology of the entire area.Geological and regional-scale literature reviews, together with acquired site-specific data, were considered as a starting point.Surface surveys were implemented by geoelectric tomography profiles and by geognostic investigations performed at various depths.A preliminary geological model allowed the planning of these geophysical and geognostic surveys, intended for the support of direct subsurface data and the elaboration of the final geological and hydrogeological models.
The recognition of the geometries and the stratigraphic relationships between the outcropping rocks and the lithological units around the town of Gimigliano, and more generally in the valleys of the Corace and Melito Rivers, were accompanied by macrostructural and meso-structural analysis to better evaluate the level of fracturing of the rock mass.

Geological Setting
The study area lies within a large sector of the Calabrian chain in the junction zone of the Sila Massif and the Catanzaro Strait (Figure 2).It is located along tectonic lines of a regional extension transversely crossing northern Calabria from the Tyrrhenian to the Ionian margin.The action of these transverse systems of regional faults is more evident in the sector just south of the study area, where the Lamezia-Catanzaro fault zone, e.g., [53][54][55], represents the master fault of a larger and more extensive system (e.g., Amantea-Gimigliano Fault).In this area, the chain is made up of Paleozoic crystalline rocks, Mesozoic metamorphic ophiolithiferous rocks, and Tertiary-Quaternary deposits (Figure 2).This group of metamorphic units, organized in numerous superimposed lithology associations, belong to the regional geological structure of Northern Calabria [54,56-58] (Figure 2).The geometrically lowest units are represented by Apennine Meso-Cenozoic carbonate sequences, overthrusted by alpine and prealpine crystalline and metamorphic rocks of low-mid-high metamorphic grade [59] (Figures 2 and 3).The group of tectonic units in the intermediate position includes mesozoic metaophiolites with metasedimentary cover and sedimentary rocks (Figure 3).The geolithological and structural analysis is the result of an original geological survey performed at a scale of 1:5000 after the photo-interpretative analysis of a significant portion of the hydrographic basin where the town of Gimigliano is located.In particular, the lithotype out-cropping in the Gimigliano area belongs to a tectonic-metamorphic unit package.From bottom to top, metamorphic rocks consisting of metapelites, metarenites, quartzites, and metalimestones are recognized in a small tectonic window and in boreholes at variable depth (ca., 60 m from surface) (Figure 3).The unit corresponds to the Frido Formation and lies on the carbonate Apennine units.
These rocks are tectonically covered by a complete ophiolitic sequence consisting of serpentinites and metabasites that overlap phyllites, metarenites, and locally thick intervals of crystalline limestones (Gimigliano-Monte Reventino Ophiolitic Unit, Figure 3 and, thanks to the interpretation of the results shown in Figure 4, see also the geological reconstruction of Figure 5).The ophiolite sequences are overthrusted by older metamorphic rocks, schists, and phyllites (Bagni Unit) and by paragneiss (Castagna Unit) (Figures 3 and 5).The reconstruction of the stratigraphic succession is complex because of polydeformed units in a ductile and then brittle environment.However, good superficial exposures of the various lithologies have allowed the identification of the main contacts between the rock associations and their state of deformation (i.e., texture, structure, fracturing, and faulting).To better investigate the geometries of the lithological units and the most important tectonic structures detected in the field survey, geoelectric tomography profiles were performed.The orientation is shown in Figure 3 and varies between mainly N-S trending and longitudinal to the slope.
The longest longitudinal profile (B-B' in Figures 3 and 4) is about 1 km for a maximum investigation depth of about 200 m.The resistivity model highlights a chaotic and mediumhigh resistivity distribution of the electro layers (Figure 4).The thickness of this "chaotic" horizon varies between 55 and 60 m, corresponding to the gravitational sliding surface.Below this interval, a resistivity distribution showing sub-vertical geometries with a lateral alternation of low-and high-resistivity rock volumes can be observed.
The interpretation of the geophysical data, together with the field geological analysis and the stratigraphic characterization of boreholes, allowed the implementation of the twodimensional geological-structural model of Figure 5.The section describes the geometries of vertical succession of the various lithologies, the antiformal fold, the major faults, the depth, and the geometry of the sliding surfaces of gravitational origin.

Hydrogeological Setting
A series of geognostic investigations were carried out to reconstruct the conceptual hydrogeological model of the study area, according to the indications provided by the geological model.Thirteen boreholes were drilled into the metamorphic aquifer at different depths and locations.During the spring season, the water table depths were measured.They are reported in Table 1 together with the geometrical characteristics of the drilled boreholes.The phreatic aquifer, characterized by low-sodium groundwater rich in calcium and magnesium bicarbonates, flows above the schists' formation, which constitutes the bottom of the groundwater body, whereas the main water circulation occurs within the serpentinites and metabasites formations (Figure 5).Its thickness varies in a range between 10 and 50 m, and the average water table depth stands at around 13 m with respect to the ground level.The aquifer is fed by superficial infiltration provided by rainfall and by a possible underground supply from a neighboring aquifer to the north.This last point is one of the key aspects of the study.The natural aquifer discharge is represented by the Corace river, a permanent water course flowing in the southern boundary of the area under investigation.
Slug tests were performed in no-dry piezometers in order to estimate the value of the hydraulic conductivity of the hard-rock Gimigliano aquifer (Table 2, Figure 3).The acquired data were interpreted using the Bouwer-Rice Method [60] because of the unconfined nature of the aquifer (Figure 6).This method is based on the mathematical model defined as follows: where h is the hydraulic head (L), r is the generic radial distance (L), and K v and K r are the vertical and horizontal components of the hydraulic conductivity (L/T), respectively.For the initial and boundary conditions, see [60].The analytical solution for the mathematical model defined by Equation ( 1), with the prescribed initial and boundary conditions [60], can be written as [61,62]: where r * w = r w √ K z /K r , H is the variation in the hydraulic head in the well with respect to the static condition (L), t is the total time from the start of the test (T), H 0 is the measurement of the initial displacement due to the slug immediately after the start of the test (L), b is the length of the well screen (L), r c is the radius of the well casing (L), R e is the effective radius parameter (L), r w is the effective radius of the well-filter section (L), and K z is the vertical component of the hydraulic conductivity (L/T).
If the slug test is performed in boreholes screened across the water table, Bouwer and Rice [61] recommend using an effective casing radius r c given by the following formula: where r nc is the nominal radius of the well screen (L) and n is the drainable porosity of the filter pack (dimensionless).
The springs falling in the study area were surveyed by cartographic analysis, field campaigns, and published studies.IGM cartography and geological maps of the Calabria region were used [63] along with technical cartography [64].Among all the emergencies identified at the cartographic level, those with greater importance (more conspicuous water flows provided with continuity throughout the year) were sought in the area as relevant for the research.
Water sources were mostly distributed in the central sector of the slope and landslide area, and according to their conspicuous frequency, the presence of a phreatic aquifer can be confirmed (Figures 7 and 8).The following springs were considered in the study as they are perennials and abundant (Table 3).The aquifer under investigation is mainly located in metabasite, serpentinite, and schist formations.The sources listed in Table 3, some of which have seasonal flowrate variations, are fed by the groundwater flowing in the aforementioned lithologies.Most of them, on the other hand, together with drainage galleries, have a quite constant flowrate during the year.
Two major draining systems were identified during the field campaigns consisting of underground tunnels.
The first one was built around the middle of the last century (in metabasites blocks inside the landslide body).It has a linear envelope with some small branches, likely built to intercept the aquifer.The tunnel is built in reinforced concrete, with the terminal parts (lateral branches and final branch) in stone.The intercepted flow is collected at the beginning of the tunnel toward the Corace river.The tunnel, for safety reasons, is not directly accessible.Finally, an important drainage element, even if made for another purpose, is the tunnel of the railway line (below the downtown), which crosses the landslide slope below the groundwater level, allowing consistent drainage and therefore not negligible for the global water balance.It was verified that the water is conveyed into the Fosso Scavone stream, often dry upstream, with a flowrate in the order of 10 L/s.

Mathematical Modeling
The mathematical model used to simulate the flow in the Gimigliano aquifer is based on the groundwater flow equation for a three-dimensional case [65]: where K is the hydraulic conductivity tensor (LT −1 ), h is the hydraulic head (L), n is the porosity (-), and Q represents the incoming/outgoing flow rate per unit volume of the aquifer (T −1 ).The universal groundwater model MODFLOW-2005 [66] was used to simulate flow phenomena into the aquifer.
The lateral recharge coming from the north domain boundary was estimated by implementing an iterative optimization procedure working in the following way: once all the model parameters were set, a first attempt value of the lateral recharge was defined on the boundary where the recharge is assumed to exist.The iterative optimization procedure was led by the Levenberg-Marquardt algorithm, used to solve non-linear least squares problems, and implemented with the model-independent parameter estimation and uncertainty analysis UCODE-2005 [67].This algorithm takes full control of the model, running it as many times as necessary to optimize the values of the desired parameter.Thus, the procedure modifies the initial value of the flow rate, assigned to the boundary, until the summation of deviation between calculated and observed hydraulic heads is minimized.

Hydrogeological Model
A three-dimensional heterogeneous model of the phreatic aquifer flowing into the study area, and working in steady state conditions, was developed.The model domain covers an area of about 1.5 km 2 and is based on the morphological and geological reconstruction of the hillside affected by the landslide (Figures 3 and 5).The 3-D geological model was compiled by integrating the surface geological map, with recently acquired and reconstructed vertical cross-sections (Section 2.1).These latter were used to recreate the layout of the aquifer bottom via an ordinary kriging approach.The ground level profile, instead, was realized using DEM data.The whole domain was horizontally discretized by a mesh of 25 × 25 m 2 square elements.The vertical extent of the model is covered by a single layer, it has a variable thickness (in a range between 10 to 50 m) that follows the trend of the impermeable formation constituting the bottom of the superficial aquifer, and it is designed to replicate the heterogeneity created by the permeable geological formations across the area (Figure 7).The horizontal hydraulic conductivities of the geological formations (Figures 3 and 5) are derived from Table 2.The vertical hydraulic conductivities were set to be a tenth of the horizontal ones, as has been done in many modeling applications, e.g., [68][69][70][71].Faults, and their influence on groundwater flow, were been considered in the model, as their occurrence in the circumscribed model domain is minor, and furthermore, detailed information on their characteristics is not available.
No flow conditions were assigned to the east and west boundaries, while an inlet flow condition was set at the north boundary.The MODFLOW River Package was then used to simulate the surface water/groundwater interaction at the south boundary (Figure 8).The permanent water course, involved as a boundary condition in the modeling process, is the Corace river.Its monthly average flowrate ranges between a maximum of 9.71 m 3 /s, in the winter season, and a minimum of 0.446 m 3 /s during the summertime [72].So, the variable flowrate will be proportional to the difference between the river stage and the hydraulic head in the aquifer, i.e., groundwater can leave the aquifer through the river boundary when the hydraulic head in the cell is higher than the river stage; otherwise, surface water can enter the aquifer through the river boundary when the head in the cell is below the river stage.
A natural aquifer recharge from rainfall N = 41 mm/y was homogeneously distributed on the model top according to the results of the water balance for the area under investigation.
The water balance can be described using the formula: where P i is the precipitation (L), E t i is the evapotranspiration (L), N i is the infiltration toward the aquifer (L), R i is the runoff (L), and ∆S i is the variation in the water storage of the soil (L).All the above-described quantities refer to the i-th month.The evapotranspiration E t i can be calculated with the method of Thornthwaite [73] based on the relationship existing between potential evapotranspiration and the average monthly temperature.The runoff R i was obtained using the infiltration curve method based on the Curve Number Method [74].
The water capacity of the ground, or field capacity, can be figuratively represented by a water-soaked sponge retaining its water content as long as the forces at play prevent its runoff.The water amount that does not runoff or evaporate seeps into the ground.If the ground is dry, the first part of this water is used to replenish the water reserve.This water asset will then be used by plants during the periods in which rain is not able to meet their water demand.The water reserve grows, until reaching the field capacity, during the rainy season while yielding the stored water when the net rain (P i − R i ) is less than evapotranspiration.When the water availability exceeds the potential evapotranspiration, and the water reserve reaches the maximum field capacity, the aquifer is recharged.Under these conditions, the variation in the water storage ∆S i is zero, and the value of the net vertical infiltration comes from the following relation: Drainage systems and springs, described in Section 2.2 and falling within the domain, were included into the model.Monitored wells and piezometers were also placed within the model and taken into consideration during the calibration procedure initially set-up to estimate the flowrate entering the system from the north inflow boundary condition (Figure 8).
As already said in the previous section, the automatic iterative optimization procedure was driven by UCODE-2005 [68].The hydraulic head values, calculated at the model observation points and representing the monitoring wells and piezometers, were compared with the measured data.The procedure changed the values of the flowrate until the deviation between calculated and observed values reached a minimum value.

Results and Discussion
The first modeling stage, represented by the inversion of the hydraulic head data, was performed to quantify the possible inlet of a flowrate from the northern boundary.The outcomes of the model at the end of the iterative optimization procedure depict a situation in which the calculated hydraulic head distribution always overcomes the observed one by several meters.The sole contribution provided by the rain is sufficient to create conditions for which the entire domain is completely flooded.This unrealistic scenario can only be attributed to the adopted hydraulic conductivity values, which are too low to allow the model to manage even only those forcings of known entity (i.e., rain recharge).
For this reason, the optimization procedure, described in the previous section, was run again to evaluate both the flow rate entering the northern boundary and the value of an equivalent hydraulic conductivity for the aquifer capable of ensuring the compliance with the observed hydraulic heads and even an incoming flow rate at least comparable to that measured at the springs.In this second modeling phase, the aquifer was then turned into a homogeneous medium, without changing the geometrical layout, as the number of available observed data samples was not sufficient to guide and constrain the estimation of the parameters of the previous heterogeneous system.The estimated hydraulic conductivity of the equivalent homogeneous system resulted in a value of 2.7 × 10 −6 m/s, which was two orders of magnitude higher than the previous permeabilities.This result, obtained for the homogeneous configuration, indirectly took also into account the effect of faults in the observed hydraulic head distribution adopted in the inversion process.
The reason of this difference between the slug tests' estimation and the numerical model inversion can be ascribed to different observation scales of the two estimates, local for the first and regional for the second.Metamorphic aquifers, in fact, are characterized by a network of secondary connections of diverse orders and degrees that, among the main fractures, realizes a continuous variation in the hydrodynamic properties at different scales defining the modulus and direction of the groundwater flow.Moreover, the slug tests were all carried out in the landslide body resulting in close permeability values for all of them.
The estimated equivalent hydraulic conductivity value (2.7 × 10 −6 m/s) enables a flow rate of 39.7 L/s entering the Northern boundary.This entry flow rate contemporarily fits the hydraulic heads of the monitored wells and piezometers and provides a groundwater flow even after the springs' feed (Figure 9a).In fact, the estimated entering flow rate is of the same order of magnitude as that of the measured springs outlets (Section 2.3.1).The goodness of the fit is shown in the scatterplot having a determination coefficient R 2 equal to 0.95 (Figure 9b).

Conclusions
In this work, an extensive campaign of field (hydro-geo-physical) measurements and a first modeling attempt of the Gimigliano aquifer are presented.The mathematical model was able to estimate the equivalent permeability of the Gimigliano aquifer and the presence of a lateral recharge from a deep confining aquifer.The adopted inverse approach provided a hydraulic conductivity of 2.7 × 10 −6 m/s and a flow rate entering the northern boundary of 39.7 L/s.
The estimated hydraulic conductivity was almost one order of magnitude higher than the permeabilities estimated by means of slug tests.The reason for this discrepancy between the slug tests' estimation and the numerical model inversion can be found in the different observation scales of the two estimates, very small for the first and regional for the latter.In fact, aquifers made up of metamorphic rocks include, among the main fractures, a network of secondary connections of different order and degree that determines, at different observation scales, a variation in the aquifer hydraulic properties defining the modulus and direction of the groundwater flow.Moreover, as a validation of our findings, the estimated lateral flow rate feeding the aquifer with a contribution of 39.7 L/s was of the same order of magnitude as the measured springs outlets.
So, as a conclusion, we argue that the hydraulic conductivity of metamorphic aquifers cannot be estimated only by field tests having a local investigation extent, but an inverse modeling approach is required, at a regional scale, which also takes into account the hydrological balance of the aquifer.

Figure 2 .
Figure 2. Schematic geological map of the northern Calabria.The black square in the southeastern sector indicates the location of the study area.Modified after Dijk et al. ii 2000 [53] and Tansi et al. ii 2007 [54].

Figure 3 .
Figure 3. Detailed geological map of the Corace River-Gimigliano village area.A-A' is the geological cross-section trace.In the central sector (white color), the landslide and slope deposits are mapped.

Figure 4 .
Figure 4. Model of resistivity of the B-B' (ERT) longitudinal section.

Figure 5 .
Figure 5. Geological cross-section along the longitudinal A-A' trace (Figure 3).The 2D geologicalstructural model is the result of the field survey, borehole stratigraphies, and ERT interpretation.

Figure 7 .
Figure 7. (a) Surface geological formations included into the model domain depicting zones with different hydrodynamic parameters; (b) 3D view of the model; (c) cross-section A-A' and (d) crosssection B-B' depicting the thickness layout of the model.

Figure 8 .
Figure 8. Boundary conditions together with relevant elements included in the model.

Figure 9 .
Figure 9. (a) Hydraulic head distribution at the end of the optimization procedure; (b) scattergraph of the observed and calculated hydraulic head values; (c) hydraulic distribution in cross sections A-A' and (d) B-B' also displayed in Figure 7.

Table 1 .
Geometrical characteristics of the boreholes drilled into the aquifer and water table depths measured during the spring season.

Table 2 .
Summary of the Slug Tests' results carried out in different piezometers.The table shows the lithology and the hydraulic conductivity obtained with the Bouwer-Rice method.

Table 3 .
Springs/Sources (S) detected in the study area.The reported values were observed in February 2011.