Monitoring and Modelling Interactions between the Montagna dei Fiori Aquifer and the Castellano Stream (Central Apennines, Italy)

: Groundwater is the most used water resource around the world, but due to population growth and climate change the alluvial lowland aquifers are often polluted and over-exploited. Thus, more and more frequently water managers need to shift their attention to mountain regions to identify groundwater resources for drinking purposes. This study presents a monitoring and modelling approach that allowed to quantify the inﬂow from the “Montagna dei Fiori” fractured aquifer to the Castellano stream. Continuous monitoring of ﬂow discharge and temperature during an entire hydrological year (2018–2019) at two monitoring stations along the stream allowed to discriminate between the baseﬂow (on average, 0.891 m 3 / s) and the run-o ﬀ (on average, 0.148 m 3 / s) components. A hydrogeological basin-wide numerical ﬂow model (using MODFLOW-2005) was set up using information from hydrogeological and geomechanical surveys. The model was calibrated using the daily baseﬂow observations made in the Castellano stream (R 2 = 0.75). The calibrated model allowed to quantify groundwater / surface water interactions. After an automated sensitivity analysis (using MODFLOW-2000), the recharge was found to be the most uncertain parameter, followed by the hydraulic conductivity zonation. This methodology could be applied in other mountain regions where groundwater monitoring networks are usually lacking to improve water resources management.


Introduction
Continuous monitoring of water fluxes and temperatures in stream reaches fed by groundwater are often key points to understanding water flux exchange from groundwater/surface water interactions [1,2]. Stream flow and temperature continuous monitoring are often a key point to understand water fluxes exchanges from and to aquifers in lowland river catchments [3][4][5], or by using remote sensing where monitoring well networks are scarce or absent [6]. As mentioned above, in recent literature there are many examples of surface waters/groundwater exchange fluxes determined by for example streambed temperature data, piezometric head data and simple 1D models or analytical solutions for flux calculation. Nevertheless, often only small-scale and site-specific information can be gained on the relationship between streams and aquifers without employing complex three-dimensional numerical models [7]. Although numerous models have been successfully developed for use in porous aquifers [8,9], their application in fractured aquifers is more problematic. Fractured media are generally highly heterogeneous since they are dominated by secondary porosity (along fractures apertures) and often are characterized by hierarchical permeability structures [10]. In these media, turbulent flow components may develop limiting the application of numerical simulations, since most numerical models are based on the Darcy's law, which is valid only for laminar flow conditions. To avoid such problems, the "black box" model approach has often been employed in the recent past to simulate fluxes in fractured aquifers [11]. The latter is a lumped parameter model, often referred to as surrogate model, where the spatial dimension is not accounted for. These surrogate models showed often good agreement between measured and simulated spring and/or stream discharges and they have been recently used to predict spring discharge [12]. The main advantage of such surrogate models is that data prerequisites are minimal and simulation runs are extremely fast, although no information is provided on internal fluxes and exchanges [13]. On the contrary, numerical groundwater flow models have been confirmed as indispensable tools to correctly design integrated water resources management plans in fractured aquifers [14,15]. These models can be employed to validate and refine previously conceived conceptual models, or to predict the impacts of climate change or anthropogenic activities on water resources availability. The main limitation of numerical flow models in fractured aquifers is due to the complex geometry of the fractures [16,17]. In many cases, the main assumption is that at large scales, fractured media can behave like an equivalent porous medium [18][19][20][21]. However, the assumptions of the equivalent porous medium and of laminar flow, in fractured systems characterized by preferential flow paths, are not always valid especially at small scales of investigation. Under these conditions, the discrete fractures approach is generally recommended [22], but its implementation is somehow restricted because it entails comprehensive data on the geometry development of fractured systems and their hydraulic properties. In any case, if the investigation is undertaken at the whole hydrogeological basin, as in the presented case, the equivalent porous medium approach can be still considered valid [23].
From the brief literature review above, it is clear that the choice of the model mainly depends on data availability, but also by a modelling code's assumptions and the objectives of the study [24]. This paper presents a methodological approach to quantify groundwater/surface water interactions using information from hydrogeological and geomechanical surveys, borehole logs, pumping tests and continuous monitoring of flow discharge and temperature. The main contribution of this research is the use of continuous discharge and temperature measurements to numerically quantify aquifer contribution to stream flow in a fractured system. The methodology was tested on a mountainous fractured aquifer, hosted by the Montagna dei Fiori, where the monitoring wells network is nearly absent. Nevertheless, with a surface monitoring system placed along the Castellano stream, which was suspected to drain the fractured aquifer, it has been possible to calibrate a complex numerical flow model employing the equivalent porous medium approach.

Materials and Methods
The MODFLOW-2005 numerical code [25] was used to create a numerical flow model of the whole hydrogeological basin of the "Montagna dei Fiori", which is a fractured aquifer that covers an area of about 50 km 2 . Pre-and post-processing were done using the graphical user interface Processing Modflow 8.0 [26]. In this aquifer the parameterization of the hydraulic conductivity (K) field was field was based on geomechanical surveys of the fracture network that were used to estimate the order of magnitude of K for the different hydrogeological units. Pumping tests conducted in wells located near the Castellano stream, provided detailed information on the K distribution in the proximity to groundwater/surface water interaction zone. Given the lack of complete information on the K field in the simulated domain, it was necessary to solve the problem with the equivalent porous media method, implicitly simulating the network of fractures that allows the groundwater flow.
Therefore, the input data required for the numerical model set up were: topography, geomechanical survey of fractures distribution, aperture and direction, hydrogeological survey, pumping tests in three wells, daily precipitation and continuous flow measurements on the stretch of the Castellano stream that Water 2020, 12, 973 3 of 14 drains the "Montagna dei Fiori" fractured aquifer. To correctly quantify the groundwater/surface water interaction at the catchment scale, the entire hydrostructure of the "Montagna dei Fiori" was considered.
The main steps of the implemented methodology to monitor and model the "Montagna dei Fiori" baseflow contribution to the investigated Castellano river reach are briefly described below.

Hydrogeological and Geomechanical Setting of the Study Area
To avoid the use of unknown boundary conditions near the groundwater/surface water interaction zone, the entire hydrogeological basin was studied ( Figure 1). Moreover, it is always advisable to use the real hydrogeological boundaries where possible, to better characterize changes in the aquifer outflow lines due to stresses such as pumping or periods with different recharge rates [13]. The impermeable boundaries of the hydrogeological basin of the hydrostructure of the "Montagna dei Fiori" are represented by folded marls, broken by the Apennine orogenetic activity. They cover the prevalently calcareous Triassic-Eocene formations, which are permeable by cracks and fractures. In contrast, the Oligo-Miocene formations are almost impermeable and outcrop in the south and in the center of the study area. In order to better detail the hydrogeological properties in the study area, a detailed map was produced for the Casteltrosino area ( Figure 2) associating the lithologies with the hydraulic properties investigated in the large scale survey, to accurately localize the compressive faults and thrusts, the pumping wells, the Castellano stream and an artificial lake formed by a hydroelectric dam.
The conceptual model consists of a confined/unconfined aquifer complex with a fractured calcareous matrix that is layered with hydrogeological units made up of finer materials (marls). This feature creates a series of semi-confined aquifers that become a single unconfined aquifer only near the tectonic window, located at the center of the ridge of the Montagna dei Fiori (Figure 1). The thickness of the non-consolidated alluvial sediments is relevant only within the artificial lake "Lago di Casette" and decreases near the rocky slopes; the bedrock of the aquifer system is given by the Burano anhydrites located on average, at −1200 m ASL.
Six stations were chosen for the survey of the geomechanical data in the area, to achieve a rough estimate of the equivalent permeability through the analysis of joints, fissures, fractures and any other discontinuity found in the rocks, and to have a general overview of the cracking condition in the rocks near the anticline axis. The geomechanical data were used to obtain the Schmidt stereoplots (not shown), indicating layouts and features of the discontinuities, from which the values of the maximum hydraulic conductivity (Kmax) and the equivalent hydraulic conductivity (Keq) were derived for each hydrostratigraphic unit, applying the Snow method [27]. Kmax and Keq were particularly high (Table 1) indicating a wide opening of the fractures due to the tensional release near the outcropping zones. This influences the infiltration in the unsaturated zone, focusing the majority of the recharge near the hinge of the anticline where the cracking system is wide open, and is thus capable of transferring a high volume of infiltrated water towards the saturated part of the aquifer [28]. Here, the groundwater flow takes place in fault and fracture systems, mostly oriented in Apennine direction. A geomechanical analysis and hydrogeological survey allowed the identification of the hydrogeological boundaries: the overthrust plain to the east and the normal extensive fault (160 • N, western dipping) to the west.
dei Fiori" basin, other than those recently implemented by the CIIP Water Company (Figure 2) drilled at mean depths of 300 m for drinking water supplies in case of extreme droughts. This is the only core log information available for the whole hydrogeological basin (not shown).

Numerical Model Set-Up
The translation of the conceptual model into the numerical model occurred through successive refining stages. The starting point was a simplified reconstruction of the boundary conditions and stresses, namely recharge rate and aquifer-streams interconnection. The second stage was the reconstruction of the three-dimensional distribution of K values via the combined information gained by the pumping tests performed by CIIP in 2018 and analyzed with the Cooper-Jacob method (mean K value of 1.8 × 10 −3 m/s) and the geomechanical surveys (Table 1), which led to the realization of the steady-state model. Furthermore, some literature values were known for the above-mentioned hydrogeological units [32,33].
The domain of the model had an extension of 14,000 m × 16,000 m and had been divided into 140 columns and 160 rows. The active cells area, in which groundwater flow and heads were computed, was of 47.91 km 2 ; while the active cells volume was of 60,968 × 10 6 m 3 with a maximum model thickness of approximately 2000 m and a maximum model depth of −600 m ASL. The vertical discretization was completed through the subdivision into 4 layers of variable thickness, following the lithostratigraphic reconstruction ( Figure 3). In Figure 3 shows the final calibrated K field used in both the steady state and transient model.
Five hydrostratigraphic units were identified based on the hydrogeological survey: (i) the Scaglia calcareous complex, which outcrops in most parts of the study area, (ii) the Marne a Fucoidi marl complex, that acts as aquitard in fault zones characterized by relaxing kinematics, (iii) the Maiolica calcareous complex characterized by a good hydraulic conductivity due to fracturing, (iv) the Jurassic siliceous calcareous complex which acts as an aquiclude and (v) the Massiccio calcareous complex characterized by a high hydraulic conductivity due to fracturing.
In this area, the recharge rate is particularly elevated since the soils are permeable, although the high terrain slope can generate elevated surface runoff. Regional literature studies put the upper limit of the recharge rate at about 600 mm/y [29,30]. No existing wells have been mapped in the "Montagna dei Fiori" basin, other than those recently implemented by the CIIP Water Company ( Figure 2) drilled at mean depths of 300 m for drinking water supplies in case of extreme droughts. This is the only core log information available for the whole hydrogeological basin (not shown).

Monitoring Strategy
To obtain hydrological data on the investigated basin, a network of rain-gauge stations and hydrometric stations was installed. These data were analyzed in conjunction with the data derived from the existing monitoring network of the Civil Protection Agency of Marche Regional Authority ( Figure 1).
Rain-gauge stations were positioned at different elevations along an ideal transect, from about 280 to 1430 m ASL. Hydrometric stations were positioned at the aquifer boundaries (upstream, where the discharge begun to increase, and downstream; Figure 2), in order to obtain the contribution of the hydrostructure to the Castellano stream. The level and temperature sensor probes (TD Diver and Barologger, Eijkelkamp Soil & Water), were placed on the stream bed within a perforated steel tube (2.5 cm internal diameter) reaching the middle of the stream to be sure that they were always submerged with water. Moreover, every 3 months from 2017 to 2019, manual flow measurements were performed using a flow meter (OTT C31, OTT Hydromet) integrating the flow velocities along the monitored stream section.
The average daily flow rates of the upstream section were subtracted from those of the downstream section to obtain the flow rate variation in the studied stretch. The obtained flow rates were further processed using the base flow index (BFI) software [31] to determine the contribution of the aquifer and to exclude the surface runoff component.

Numerical Model Set-Up
The translation of the conceptual model into the numerical model occurred through successive refining stages. The starting point was a simplified reconstruction of the boundary conditions and stresses, namely recharge rate and aquifer-streams interconnection. The second stage was the reconstruction of the three-dimensional distribution of K values via the combined information gained by the pumping tests performed by CIIP in 2018 and analyzed with the Cooper-Jacob method (mean K value of 1.8 × 10 −3 m/s) and the geomechanical surveys (Table 1), which led to the realization of the steady-state model. Furthermore, some literature values were known for the above-mentioned hydrogeological units [32,33].
The domain of the model had an extension of 14,000 m × 16,000 m and had been divided into 140 columns and 160 rows. The active cells area, in which groundwater flow and heads were computed, was of 47.91 km 2 ; while the active cells volume was of 60,968 × 10 6 m 3 with a maximum model thickness of approximately 2000 m and a maximum model depth of −600 m ASL. The vertical discretization was completed through the subdivision into 4 layers of variable thickness, following the lithostratigraphic reconstruction ( Figure 3). In Figure 3 shows the final calibrated K field used in both the steady state and transient model.
The top of the first layer was not horizontal, but it was mimicking the real topography, by re-interpolating a 100 × 100 m resolution digital terrain model (DTM) on the MODFLOW-2005 grid. Where possible, the distance between the calculation node of the cells belonging to a given layer and the calculation nodes of the cells belonging to the adjacent layers was set not greater than 50% of the thickness of the layer, to facilitate the numerical model convergence. The base of the model domain (bottom of layer 4) was interpolated using the hydrogeological sections produced during the hydrogeological survey. The constant head package was used to simulate the piezometric head in the Casette Lake ( Figure 2); this lake and the Castellano stream were considered the hydrogeological basin outflow, while all the model cells outside the hydrogeological basin ( Figure 1) were set as no flow boundary. The top of the first layer was not horizontal, but it was mimicking the real topography, by reinterpolating a 100 × 100 m resolution digital terrain model (DTM) on the MODFLOW-2005 grid. Where possible, the distance between the calculation node of the cells belonging to a given layer and the calculation nodes of the cells belonging to the adjacent layers was set not greater than 50% of the thickness of the layer, to facilitate the numerical model convergence. The base of the model domain (bottom of layer 4) was interpolated using the hydrogeological sections produced during the hydrogeological survey. The constant head package was used to simulate the piezometric head in the Casette Lake ( Figure 2); this lake and the Castellano stream were considered the hydrogeological basin outflow, while all the model cells outside the hydrogeological basin ( Figure 1) were set as no flow boundary.
The Castellano stream and the Salinello stream, which represent the only perennial streams in the hydrogeological basin, were simulated using the river package. The input values entered were the riverbed mean K value that was estimated by the Hazen formula [34] (equals to 1.0 × 10 −4 m/s) and the height of the riverbed level which was obtained from the DTM and the field survey. Clearly, the mean K value is an assumption that does not take into consideration the various degree of heterogeneity that are likely to be present along the stream reach [8], but this was the best estimate that was available. The river package made use of cells in which the hydrometric river stage was specified by the user (here two stages for the Castellano stream and three for the Salinello stream at monitored points); a linear interpolation was applied between the river cells with user specified values. The horizontal flow barrier (HFB) package was used to simulate the compressive faults, thus including large anisotropy in the K field at the basin scale.
The transient simulation was run for the period in which it was possible to continuously monitor the flow of the Castellano stream (October 2017-March 2019). The time was divided into ten stress periods to reproduce the main recharge events; the length of the periods was inferred from the peak analysis recorded in the two continuous measurement hydrometric stations along the Castellano stream. Furthermore, the average hydrometric levels of the Castellano stream in each of the ten stress periods was derived from the continuous monitoring data of the two hydrometric stations. An automated calibration process was carried out using the numerical inverse code, MODFLOW-2000 (parameter estimation and sensitivity analysis) [35]. This process was performed on the following parameters: horizontal K of the four hydrogeologic units, hydraulic conductance of riverbed and recharge. Prior information on the isotropic distribution of Kx and Ky values and on the ratio between horizontal and vertical K values were used in order to constrain parameters during estimation; the The Castellano stream and the Salinello stream, which represent the only perennial streams in the hydrogeological basin, were simulated using the river package. The input values entered were the riverbed mean K value that was estimated by the Hazen formula [34] (equals to 1.0 × 10 −4 m/s) and the height of the riverbed level which was obtained from the DTM and the field survey. Clearly, the mean K value is an assumption that does not take into consideration the various degree of heterogeneity that are likely to be present along the stream reach [8], but this was the best estimate that was available. The river package made use of cells in which the hydrometric river stage was specified by the user (here two stages for the Castellano stream and three for the Salinello stream at monitored points); a linear interpolation was applied between the river cells with user specified values. The horizontal flow barrier (HFB) package was used to simulate the compressive faults, thus including large anisotropy in the K field at the basin scale.
The transient simulation was run for the period in which it was possible to continuously monitor the flow of the Castellano stream (October 2017-March 2019). The time was divided into ten stress periods to reproduce the main recharge events; the length of the periods was inferred from the peak analysis recorded in the two continuous measurement hydrometric stations along the Castellano stream. Furthermore, the average hydrometric levels of the Castellano stream in each of the ten stress periods was derived from the continuous monitoring data of the two hydrometric stations. An automated calibration process was carried out using the numerical inverse code, MODFLOW-2000 (parameter estimation and sensitivity analysis) [35]. This process was performed on the following parameters: horizontal K of the four hydrogeologic units, hydraulic conductance of riverbed and recharge. Prior information on the isotropic distribution of Kx and Ky values and on the ratio between horizontal and vertical K values were used in order to constrain parameters during estimation; the objective function consisted of all the available flow observations. R 2 , the Nash-Sutcliffe model efficiency coefficient (NSE) and the absolute mean error (MAE) were used to infer if model calibration was satisfactory.

Monitoring Results
The observed inflow and outflow from the Castellano stretch that directly intercepts the Montagna dei Fiori aquifer shows a marked increment of the stream flow in this small portion of about 1 km length Figure 4, the increase in stream discharge due to aquifer recharge reached a maximum of about 1.2 m 3 /s in the first hydrological year and 0.7 m 3 /s in the second year. During the monitoring period, several flood events occurred, always in concomitance with the precipitation falling within the basin. In the left part of the graph in Figure 4, the increase of the stream discharge due to the aquifer recharge is particularly evident (see Castellano Downstream red line, from November 2017 to mid-April 2018). objective function consisted of all the available flow observations. R 2 , the Nash-Sutcliffe model efficiency coefficient (NSE) and the absolute mean error (MAE) were used to infer if model calibration was satisfactory.

Monitoring Results
The observed inflow and outflow from the Castellano stretch that directly intercepts the Montagna dei Fiori aquifer shows a marked increment of the stream flow in this small portion of about 1 km length Figure 4, the increase in stream discharge due to aquifer recharge reached a maximum of about 1.2 m 3 /s in the first hydrological year and 0.7 m 3 /s in the second year. During the monitoring period, several flood events occurred, always in concomitance with the precipitation falling within the basin. In the left part of the graph in Figure 4, the increase of the stream discharge due to the aquifer recharge is particularly evident (see Castellano Downstream red line, from November 2017 to mid-April 2018). The temperature behaviour is even more explicit than the observed stream flow, in fact while the temperature in the Castellano Upstream Station was clearly following the daily atmospheric fluctuations, the temperature in the downstream monitoring station of the Castellano stream did not vary except for four storm events when a larger component of the runoff contributed to a decrease the stream water temperature ( Figure 5). This is a clear indication that the aquifer contribution in the monitored Castellano stretch is extremely significant as pointed out in other environmental settings by other studies [36,37]. The temperature behaviour is even more explicit than the observed stream flow, in fact while the temperature in the Castellano Upstream Station was clearly following the daily atmospheric fluctuations, the temperature in the downstream monitoring station of the Castellano stream did not vary except for four storm events when a larger component of the runoff contributed to a decrease the stream water temperature ( Figure 5). This is a clear indication that the aquifer contribution in the monitored Castellano stretch is extremely significant as pointed out in other environmental settings by other studies [36,37]. Figure 6 shows the baseflow calculated via the BFI model using as input value the difference between the upstream and downstream observed flows in the Castellano stream. BFI eliminates the run-off component (flow peaks) using a simplified algorithm. This approach can be considered valid in this simple case, since the Castellano stream is the only stream reach of this area, but for complex cases like in the presence of a dense stream network or in urban environments, a distributed watershed model should be employed to accurately calculate the baseflow [38]. The baseflow feeding the Castellano stream over the monitored period is at least 400 l/s with large variations induced by rainfall events, see for example the 120 mm event of November 2017 that triggered a large baseflow increase. For instance, from November 2017 to mid-April 2018 around 1500 l/s were discharged by the aquifer into the Castellano stream. Accordingly, this period of high baseflow coincides with a particularly rainy meteorological period, in fact as many as 580 mm were recorded by the meteorological station of the Civil Protection of the Marche Region, while in the same period of the following year only 314 mm were recorded by the same station.   Figure 2) and the calculated BFI stream baseflow. Rainfall for the same period is also shown. Figure 6 shows the baseflow calculated via the BFI model using as input value the difference between the upstream and downstream observed flows in the Castellano stream. BFI eliminates the run-off component (flow peaks) using a simplified algorithm. This approach can be considered valid in this simple case, since the Castellano stream is the only stream reach of this area, but for complex cases like in the presence of a dense stream network or in urban environments, a distributed watershed model should be employed to accurately calculate the baseflow [38]. The baseflow feeding the Castellano stream over the monitored period is at least 400 l/s with large variations induced by rainfall events, see for example the 120 mm event of November 2017 that triggered a large baseflow increase. For instance, from November 2017 to mid-April 2018 around 1500 l/s were discharged by   Figure 2) and the calculated BFI stream baseflow. Rainfall for the same period is also shown. Figure 6 shows the baseflow calculated via the BFI model using as input value the difference between the upstream and downstream observed flows in the Castellano stream. BFI eliminates the run-off component (flow peaks) using a simplified algorithm. This approach can be considered valid in this simple case, since the Castellano stream is the only stream reach of this area, but for complex cases like in the presence of a dense stream network or in urban environments, a distributed watershed model should be employed to accurately calculate the baseflow [38]. The baseflow feeding the Castellano stream over the monitored period is at least 400 l/s with large variations induced by rainfall events, see for example the 120 mm event of November 2017 that triggered a large baseflow increase. For instance, from November 2017 to mid-April 2018 around 1500 l/s were discharged by Figure 6. Flow rate difference between the two instrumented sections of the Castellano stretch (upstream and downstream hydrographic stations in Figure 2) and the calculated BFI stream baseflow. Rainfall for the same period is also shown.

Steady State Model Results
During the model calibration process, the initial estimates of K values for the different hydrogeological units were changed within the limits of physical consistency, to obtain calculated piezometric heads values (approximately 312 m ASL for all the wells) comparable with the ones measured (approximately 310 m ASL for all the wells) in the wells drilled by CIIP in 2018 near Casteltrosino (see Figure 2 for locations). These were the only piezometric head data available to constrain the steady state model.
The effective rainfall (recharge) was simulated with the recharge package, the estimate of the effective annual infiltration was initially set to an average value of 600 mm/y (30% of precipitation), but the results showed a great discrepancy with the observed flow measurements in the Castellano stream. In fact, they were orders of magnitude higher (8000-10,000 l/s) than those estimated by discontinuous flow rate measurements (500-900 l/s). The recharge value was then gradually reduced, reaching a mean value of about 378 mm/y (19% of precipitation). The water table in the steady state model had a variable depth with respect to the ground level, very deep near the topographic peaks and approximately above ground level near the Castellano stream (Figure 7). The aquifer was characterized by higher hydraulic gradients in the northern portion of the flow domain near to the Castellano stream and the Casette Lake. The hydraulic head gradient on average was equal to 0.5% , therefore extremely low due to the elevated K values of the modeled hydrogeological units. The main groundwater flow direction was SSE-NNO along the ridge of the Montagna dei Fiori, locally influenced by variations in K and by the Salinello stream, which is located at a higher altitude than the water table, thus is only feeding and never draining the aquifer system. Furthermore, the fault system did not allow groundwater to flow through the fault planes, but instead groundwater flow was conveyed along the fault planes (Figure 7).

Transient State Model Results
The results of the transient simulation highlighted that the piezometric fluctuations induced by the different recharge rates applied during the different stress periods provoked changes in the Castellano stream baseflow. A piezometric peak was reproduced in January 2018, while the minimum occurred in February 2019 (piezometric maps not shown). However, the general trend of the piezometric heads remains, which moves sub-parallel to the main fracturing direction along the ridge of the Montagna dei Fiori, and by then it realigns along the direction of the Castellano valley.
The calibration of the model has reached a good degree of accuracy concerning the calculation of the flow rates towards the Castellano stream, which fluctuate mimicking those recorded by the continuous monitoring stations, as shown in Figure 8. Moreover, also the scatter diagram for the direct comparison between the observed and the calculated flow rates, gives acceptable R 2 (0.75), NSE The water balance is shown in Table 2. The outflow of the aquifer towards the Castellano stream is approximately 530 l/s and takes place in the central portion of the stretch here simulated, while the recharge from rivers to the aquifer is 67 l/s (calculated flow) and it is due to the Salinello stream. These data coincide with the flow measurements along the Castellano stream and the Salinello streams during the field surveys. As pointed out before, the annual recharge calculated in the study area is about 19% of the average rainfall observed over the whole hydrogeological basin (1990 mm/y) considering all the meteorological stations available, which was very similar to the mean annual precipitation for the period 1991-2020 (approximately 1800 mm/y) observed in this this area [32]. The estimated mean value of 380 mm/y is consistent with the recharge values obtained in similar hydrogeological settings [33]. The mean recharge rate multiplied by the area of the active cells in the numerical model is equal to an inflow of about 575 l/s ( Table 2).
Considering all the components of the water balance, the inflows and outflows in the study area are around 640 l/s, a result that is consistent with the water balances drawn up with simplified balance equations in similar settings [33,35].
Unfortunately, since the system is very complex from a structural point of view, with faulted hydrogeological units and boundaries that cannot be easily reconstructed in depth, it is difficult to assign consistent limits of reliability to the simulation without observed head data to calibrate the model. Thus, a transient flow model was built up and calibrated versus observed flows as reported in the next section.

Transient State Model Results
The results of the transient simulation highlighted that the piezometric fluctuations induced by the different recharge rates applied during the different stress periods provoked changes in the Castellano stream baseflow. A piezometric peak was reproduced in January 2018, while the minimum occurred in February 2019 (piezometric maps not shown). However, the general trend of the piezometric heads remains, which moves sub-parallel to the main fracturing direction along the ridge of the Montagna dei Fiori, and by then it realigns along the direction of the Castellano valley.
The calibration of the model has reached a good degree of accuracy concerning the calculation of the flow rates towards the Castellano stream, which fluctuate mimicking those recorded by the continuous monitoring stations, as shown in Figure 8. Moreover, also the scatter diagram for the direct comparison between the observed and the calculated flow rates, gives acceptable R 2 (0.75), NSE (0.53) and AME (0.146 m 3 /s), even if the model tends to overestimate the flow rates entering the Castellano stream during the high discharge periods. This is likely due to the fact that the simulated medium is fractured, while the model corresponds to an equivalent porous medium. Castellano stream during the high discharge periods. This is likely due to the fact that the simulated medium is fractured, while the model corresponds to an equivalent porous medium.

Sensitivity Analysis Results
After conducting the model sensitivity analysis to the perturbation of the input parameters with MODFLOW-2000, it was clear that the most influencing parameters were the recharge values and the conductance value of the Castellano stream ( Table 3). Given that the shape and geometry of the Castellano stream are known, this means that the K of the riverbed is the most influencing parameter after the recharge values in stress periods 2 and 3. The model sensitivity to the perturbation of the K values for the layers 1 and 2 is less pronounced and decreases further when K values set with the HFB package are also perturbed. This means that the geomechanical and hydrogeological

Sensitivity Analysis Results
After conducting the model sensitivity analysis to the perturbation of the input parameters with MODFLOW-2000, it was clear that the most influencing parameters were the recharge values and the conductance value of the Castellano stream ( Table 3). Given that the shape and geometry of the Castellano stream are known, this means that the K of the riverbed is the most influencing parameter after the recharge values in stress periods 2 and 3. The model sensitivity to the perturbation of the K values for the layers 1 and 2 is less pronounced and decreases further when K values set with the HFB package are also perturbed. This means that the geomechanical and hydrogeological characterizations were of sufficient quality to allow a unique K field reconstruction. The high recharge and river conductance sensitivities are not surprising, since these parameters drive the groundwater/surface water interactions in this numerical model; besides, they tend to be the most uncertain parameters in large scale models [39,40]. The final riverbed K value was 1.7 × 10 −4 m/s, which was similar to the initial estimate.

Conclusions
The continuous monitoring of the Castellano stream provided insight into the relationships between surface waters and the fractured aquifer. In particular, the temperature monitoring provided an independent information on the contribution of groundwater discharge to the Castellano stream. However, it must be stressed that the temperature information only provided additional qualitative information on the considerable contribution of groundwater to streamflow but was not considered for the baseflow analysis and neither used for further flux quantifications in this study. The transient flow model was able to reproduce the stream baseflow reasonably well. Nevertheless, sub-reach interactions could not be assessed with the available dataset; in fact, only interactions due to large-scale water level differences between aquifer and stream were taken into account but not the processes occurring due to local flow processes along the streambed. The automated sensitivity analysis, performed with MODFLOW-2000, showed that recharge was the most uncertain parameter, followed by hydraulic conductivity zonation and river conductance. Thus, these parameters should be better investigated in future studies in similar mountain regions where groundwater monitoring networks are usually lacking. Nevertheless, the presented monitoring and modelling approach could be employed to quantify the groundwater/surface water interactions in analogous hydrogeological settings.