Integration of photogrammetry from unmanned aerial vehicles, field measurements and discrete fracture network modeling to understand groundwater flow in remote settings: test and comparison with geochemical markers in an Alpine catchment

Fast and cost-effective techniques for hydrogeological modeling are of broad interest for water resources exploitation, especially in remote settings, where hydrogeological measurements are difficult to perform. Unmanned aerial vehicles (UAV)-based techniques are potentially useful for these aims, but their application is still limited. In this study, a field-based approach and UAV-based approach are integrated for the computation of a discrete fracture network model of a fractured aquifer in the Central Alps. Then, calculated directions of the hydraulic conductivity components were compared with a geostatistical analysis of geochemical markers from sampled spring waters, to infer a conceptual model of groundwater flow. The comparison of field-based and UAV-based fracture measurements confirmed a good matching for fracture orientations and recognized a more reliable estimation of fracture dimensions for the UAV-based dataset. Nonetheless, an important variable for hydrogeological modeling—fracture aperture—is not measurable using UAV, as this requires field measurements. The calculated directions of the main conductivities fit well with the analyzed geochemical markers, indicating the presence of two partially separated fractured aquifers and describing their possible groundwater flow paths. The adopted integrated approach confirms UAV-based measurements as a potential tool for characterization of fracture sets as the input for hydrogeological modeling and for a fast and effective surveying tool, reducing time and cost for other following measurements.


Introduction
The increase of human population has been inevitably affecting water demand worldwide, and concurrent climate change has been changing natural water abundance also in previously unaffected areas (Hilberg 2016). Quantification and management of water resources is thus becoming a critical issue in geological and environmental sciences, and increasing attention is focused on high mountain areas (e.g., Hilberg and Riepler 2016;Idrizovic et al. 2020) where fractured aquifers are particularly hard to be investigated (Berkowitz 2002;Voeckler and Allen 2012;Margat and Van der Gun 2013;Welch and Allen 2014;Hilberg and Riepler 2016).
Classical hydrogeological investigation methods such as groundwater head mapping and pumping tests, are suitable in urban, plain and coastal settings (e.g. Moench 1997;Castorina et al. 2013;Biswas et al. 2014). However, these approaches are unsuitable in remote areas where the morphology can affect well performance and durability-e.g. caused by steep and unstable slopes in mountainous settings which can damage piezometers-, the limited soil layer thickness complicates drilling, and the limited infrastructures can make field operations difficult. Also, especially in developing countries, borehole drilling results in an expensive tool for water source exploration (Vallet et al. 2015).
Indirect methods are valid alternatives to characterize groundwater flow in these settings. These methods include hydrochemistry surveys using natural and artificial tracers found in spring waters (Cervi et al. 2012;Hilberg 2016), fracture analysis at outcrop scale (Chesnaux et al. 2009;Voeckler and Allen 2012), and remote sensing techniques by means of unmanned aerial vehicles (UAV), to observe lineaments or to directly derive three-dimensional (3D) fracture measurements from virtual outcrops. The latter is of growing interest in environmental geosciences, as it is less time-consuming and more cost-effective. The fields of application are wide and they include analysis of: slope stability and landslides (Fernández et al. 2016;Salvini et al. 2017; Barlow et al. 2017;Casagli et al. 2017); structural geology (Cawood et al. 2017;Nesbit et al. 2018); cultural heritage and archeology (Spanò et al. 2018); geomorphology (Silva et al. 2017); bathymetry (Bandini et al. 2018); and reservoir characterization (Miranda et al. 2018). Moreover, UAV setups permit the application of multispectral imaging and other tools for data collection (e.g., thermal sensors). However, the application of UAV-based techniques to the investigation of fractured aquifers (DeBell et al. 2016;Bisdom et al. 2017a;Massaro et al. 2018) is still a relatively unexplored field compared to other geological applications (Giordan et al. 2020).
The main advantages of UAV-based remote sensing reside in time-saving data collection along with fractures and lineaments detection, including data on features/shape orientation, dimensioning and spatial distribution (Bisdom et al. 2017b). Nonetheless, this approach presents some important limitations such as the impossibility to measure fracture apertures (Bisdom et al. 2014;Seers and Hodgetts 2014;Welch and Allen 2014). Consequently, an integration of field-and UAV-collected data could meet the potential of both the approaches and a further comparison can be proposed with other proxy data such as geo-chemical makers (e.g., Gambillara et al. 2013;Vallet et al. 2015;Achtziger-Zupančič et al. 2017a).
In this paper, such an integrated approach is tested in a small Alpine catchment, where groundwater is hosted by a fractured rock aquifer, mostly covered by slope and glacial deposits. Structural datasets from traditional mapping and UAV-derived datasets are compared and integrated. Then, discretized tensors for rock permeability are computed by means of a discrete fracture network approach. Pros and cons of field techniques and UAV data collection techniques are compared to investigate the factors affecting the model outputs, and the conceptual model of water flow in the catchment is compared with geochemical markers of springs.

Study area
The study area is the Ventina Valley, located in the Italian Central Alps and covering ca. 3 km 2 with elevations ranging between 2,450 and 1,960 m above sea level (asl). It is characterized by a cold and temperate climate, with a mean annual temperature of 2°C and a mean annual precipitation of 1,123 mm (Isotta et al. 2014

Tectonic and geochemical framework
The northern sector of the study area is crosscut by an E-W striking subvertical fault (Pirola fault; Fig. 1). The fault separates different alpine units, which have been alternatively interpreted as fully included in the Pennidic domain with extensional allochthonous units of the Adria continental margin (Bonsignore et al. 1971;Trommsdorff and Evans 1977) or as belonging to the Pennidic Suretta nappe in the south and the Austroalpine Margna nappe in the north (Coward and Dietrich 1989;Schmid et al. 2004 (Bedogné et al. 1993;Trommsdorff et al. 2005). In the north-eastern area, the lithology consists of albitic and chloritic paragneiss (Bonsignore et al. 1971 Bedogné et al. 1993).
To the south, the dominant lithology is represented by greenschist facies serpentinites, which are metamorphic rocks derived from igneous Mg-Fe-rich peridotite. The main minerals include aggregated or sheeted antigorite [(Mg,Fe) (Bonsignore et al. 1971), as well as accessory minerals containing heavy metals such as Ni and Cr (Kierczak et al. 2008;Morrison et al. 2015;Binda et al. 2018).
In the proximity of the glacier tongue, an ophicarbonate lens is outcropping. This NW-SE striking tabular body is ca. 400 m thick and is exposed for a length of 6 km (Bonsignore et al. 1971;Trommsdorff and Evans 1977). The Ventina ophicarbonates are characteristically inhomogeneous metamorphic rocks with a brecciated texture, consisting of mmto m-sized fragments of serpentinite, embedded in a fine-to medium-grained matrix of calcitic composition (CaCO 3 ) with associated talc veins (Pozzorini and FruhGreen 1996).

Hydrology
The study area consists of two main hydrological basins: (1) Ventina glacier basin, covering ca. 11 km 2 (of planimetric surface), where an ice tongue actively supplies the Ventina River and (2) the adjacent Pirola Lake basin, recharged by atmospheric precipitation, rock glaciers and snowfields, extending over ca. 1.5 km 2 . The Pirola Lake basin is divided into two subbasins (Pirola 1 and Pirola 2 in Fig. 1), separated by a roughly NW-SE striking ridge of rouche mountonné. Pirola basin is endorheic but lake water potentially flows out through an artificial spill point at the NW tip of the lake (Radice et al. 2016).

Hydrogeology
Two main phreatic aquifers can be recognized in the study area. In the glacial and talus deposits extensively covering the bedrock, a relatively fast groundwater flow, fed by precipitation and snow melt, is mainly controlled by the amount of fines in the sediment matrix. The juxtaposition of diverse sedimentary facies determines the presence of local aquitards that locally force the groundwater to emerge (Hilberg and Riepler 2016). Several well-preserved frontal and lateral moraines are present along the Ventina Valley and are ascribable to the Little Ice Age between ca. 1,300 and 1,800 AD; (Matthews and Briffa 2005). Older lateral moraines probably date back to the Last Glacial Maximum (i.e., around 30-18 ka BP, Ivy-Ochs et al. 2008) and late-glacial stages (Trommsdorff et al. 2005). Locally, these constitute an aquitard, forcing the groundwater to rise and flow from the Pirola shallow aquifer to the Ventina Valley aquifer. Thus, in this shallow aquifer, water flow is mainly driven by the spatial distribution of glacial and periglacial morphologies and by the availability of water, mainly by snow melt.
A second aquifer is hosted by the fractured bedrock, with diffuse point springs characterized by a continuous a b Fig. 1 a Regional setting, and b geological, morphological and hydrological setting (coordinates in WGS 84 format) of the study area, including locations of: measuring sites for the field-based model, virtual outcrops used for the UAV-based model, and sampling points for water chemical analyses flow, at least during the summer season. The primary permeability of the local basement rocks is generally below 10 −16 m 2 (Fisher 1998;Zharikov et al. 2005;Achtziger-Zupančič et al. 2017b); thus, it is presumable that all the groundwater flow occurs within fractures. Figure 2 displays the workflow for the permeability tensor computation from the integrated DFN model, derived from classical structural mapping and UAV remote sensing (Lang et al. 2014) and for its qualitative comparison with the interpolated geochemical analyses.

Geological measurements and modeling
Outcrop-scale data collection Geological maps were first analyzed and 21 representative natural outcrop locations ( Fig. 1) of lithologically and structurally distinctive rock groups were selected for structural mapping (Welch and Allen 2014).  The geometrical properties of the fracture networks were measured together with the orientations. For the latter, the Android app FieldMOVE Clino by Petroleum Experts was used and the results were compared with a standard geological compass. Several measurements for each fracture cluster were collected, then features and statistics were established for the geometrical properties.

Flow modeling
The georeferenced datasets were then imported in the MOVE structural software by Petroleum Experts for structural fracture modeling and analysis, with the Fracture Analysis module. Single outliers were removed after Q-Q plot analysis and application of Fisher's statistic (Fisher 1953). Orientation data were clustered into fracture sets using a K-means clustering procedure (Tokhmechi et al. 2011). Average orientations and standard deviations were calculated for each of the fracture sets.

UAV flight, virtual-outcrop model generation and data measurements
Three flights were performed, partially covering the study area, by means of a commercial UAV (DJI Phantom 4 Pro Plus) equipped with a built-in high-resolution optical camera. Flights were performed on 15 September 2019 for ca. 20 min each. Over 900 images of three of the most representative outcrops in the Pirola basin (Figs. 1 and 3) were collected, in order to cover a zone presenting also a high number of measuring points. Georeferenced virtual outcrops were established for more detailed information, e.g. of the fracture length and width. Field measurements of these parameters are usually difficult due to the limited accessibility of these outcrops.
In order to reconstruct the 3D digital models, the digital photos were processed with a structure from motion (SfM) technique, developed in the software Agisoft Metashape Professional. Photos were lens-calibrated, and then processed following the best-practice workflow of the software, in order to obtain a point cloud density adequate for the study's purposes. Obtained HR point clouds typically consist of ca. 10 7 points with an average spacing of less than 1 cm at the ground (Fig. 3). Structural data were measured through the Compass plugin (Thiele et al. 2017) of the CloudCompare software (CloudCompare 2020). Each fracture plane was extracted by digitizing the intersection between the fracture trace and the topography. Compass provides a semi-automatic tracing of the intersection, based on changes in point color and edge detection. The extrapolated fracture plane results from the best-fit rectangular plane, passing through the support vertices of the digitized trace. Output data include rectangular plane orientation, minimum width and minimum length. Finally, fracture data were imported in MOVE and analyzed with the same statistical procedure as for the field-based measurements (see section 'Outcrop-scale data collection').

Fracture modelling from field and virtual-outcrop measurements and fractured aquifer characterization
A stochastic discrete fracture network (DFN) was generated, based on field and virtual outcrop data, in order to characterize the directions of anisotropy in the hydraulic conductivity within the rock volume. As input data, the model requires the average and uncertainty values of orientation, spacing and aperture for each of the fracture sets. As output, the modeling populates a volume with the three-dimensional (3D) eigenvectors of the permeability tensor.
The boundary conditions of the model are (1)   2017; Trinchero et al. 2017), (2) unsaturated conditions of the aquifer, since reliable measurements of the groundwater gradient and direction are lacking, (3) a depth-independent aperture value for the fractures and (4) a laminar water flow inside the fractures (Oda 1985). The study did not consider the relative orientation between fractures and the regional stress field, since it was considered that the main mechanism responsible for rock fracturing is post-glacial stress changes (e.g., regional-scale post-glacial rebound and/or valley side releasing) and regional stresses are strongly affected by near-surface reorientation due to the topographic effect. Nonetheless, surface parallel jointing exfoliation was not observed, indicating that the observed discontinuities are the final result of complex concurrent mechanisms. For the same reason, lacking reliable in-situ measurements of the maximum horizontal tectonic stress, a quantitative modeling of the stress conditions with depth and of their effect on fracture aperture was omitted from the model. Given all above, the results from the model were considered as valid for a limited depth of investigation (i.e., ca. 100 m below ground surface, bgs) where conditions similar to those ones observed at surface can be assumed.
In this study an integration of UAV-based and field-based data was applied, based on clustering of fracture orientation and on the relative spatial location of measurement sites. This approach allowed the researchers to overcome the common limitations deriving from populating a DFN based on field measurements alone (i.e., the limited reliability of fracture dimensions and aspect ratio).
The selected area for the model includes the northern part of the Ventina basin and the southern part of the Pirola basin, covering about 3 km 2 of planimetric surface (Fig. 1).
Firstly, the rock volume was divided into separate and discrete sub-volumes: & The study area was divided into sectors by creating Thiessen polygons (De Berg et al. 2008) around each field site (areas are coded as the measuring sites for sake of simplicity) and adding the Pirola fault as an additional element for the definition of the polygon boundaries. This is a common deterministic procedure for the spatial discretization, and the extent of the polygons is predetermined by the spatial distribution of data alone; & From every single Thiessen polygon, a geo-cellular volume was generated (cubic cells, 20 m spaced) using the digital elevation model as the top boundary and a horizontal lower boundary located at the average Ventina Valley elevation at 2,000 m asl (i.e., encompassing an elevation range of ca. 400 m).
Then, a DFN was generated for every subvolume. Firstly, the clusters of fractures from the field measurements were compared with those ones coming from the virtual outcrops. Comparison was performed based on the relative locations of field stations and virtual outcrops: all the data collected at measuring stations were compared with the closer virtual outcrop and clusters with similar orientations were considered as representative of the same fracture set (e.g., Voeckler and Allen 2012). A new average orientation and Fisher's K value were then recalculated (Table 1). In the case of inconsistences of cluster orientation, values of UAV or field locations only were used for the model, including the relative assumptions (Table 1). When the cluster of fractures presented both UAVbased and field-based measurements, the integration included spacing, dimension and aspect measures obtained from the virtual outcrops and the aperture observed at the field site (see section 'Outcrop scale vs UAV-based measurements: comparison of data input' for a critical comparison of UAV and field measurements). All the input data used for the DFN generation are reported in Table S1 of the electronic supplementary material (ESM1).
Then, the permeability tensor was calculated using the method of parallel plate by Oda (1985) describing the laminar flow of fluids along smooth fracture systems. More details regarding the DFN generation and permeability tensor calculation are included in Appendix 1.
Permeability ellipsoids were used for the evaluation of results in 3D and for representing planar anisotropies in the permeability directions (Lang et al. 2014). Since equivalent permeability values (k) presented in this study span several orders of magnitude, normalization is performed with respect to the local k max. This normalization applies to all presented plots and implies that ellipses and ellipsoids may be only used to compare orientation and degree of anisotropy of equivalent permeability, but not absolute magnitudes. Finally, for the comparison with geochemical markers, the 3D permeability ellipsoids were projected onto 2D maps to evaluate plane view anisotropies.

Water-samples collection and chemical analyses
The major ions and trace elements composition of water is strongly correlated with the aquifer bedrock geochemistry (Achtziger-Zupančič et al. 2017a). This is especially true in mountain catchments, where the soil is scarcely developed (Tranter 2003;Binda et al. 2018). In this study, major ions and trace elements data were processed as a comparison with the conceptual model obtained from fracture analysis.
Water samples were obtained monthly in a total of 10 sampling campaigns during summer 2014, 2015 and 2016. The samples were collected from June to October because of the high-mountain climate conditions, making sampling impossible during winter.
The points selected for chemical analyses in this study partly takes note of the points selected for trace elements source apportionment in the study area by Binda et al. 2020a. In this study, in order to observe chemical quality of groundwater only, solely springs from fractures and from talus and slope deposits were analyzed. Water derived from the direct melt of the Ventina glacier was not considered here.
More details about chemical analyses of waters are reported in Binda et al. (2020a).
Water chemical data were analyzed by principal component analysis (PCA). This multivariate statistical tool well highlighted the correlations between variables and the main geochemical trends (Gambillara et al. 2013;Giussani et al. 2016;Rosen et al. 2018).

Results and discussion
The comparison of different data-collection techniques enables the best input data to be set up for the integrated permeability model. Additionally, the chemical analysis of springs in the study area allows the observation of reliable geochemical markers of water flow, making it possible to infer a conceptual model of the aquifer in the study area.

Outcrop-scale vs UAV-based measurements: comparison of data input
In this section, the results obtained by the different data measurement techniques regarding fracture orientations and dimensions are compared, in order to elucidate the choices made for the integrated model construction. The two data input techniques, as observed in the 'Materials and methods' section, clearly require some assumptions, which can affect the final flow outputs (Bisdom et al. 2014;Seers and Hodgetts 2014;Welch and Allen 2014). The comparison of the strength and the flaws of these techniques will permit a more consistent integration of their measurements, obtaining a more representative flow model.

Fracture orientations
A good consistency is observable, combining most of the data from field survey with virtual outcrops, highlighting the potential for a fast and effective data collection based on UAVderived 3D models (Bisdom et al. 2017a;Akara et al. 2020), even spanning over different scales (Giuffrida et al. 2020). Figure 4 shows a comparison between data derived from virtual outcrops and from field measurements. Four subareas were identified and three of them were associated to virtual outcrops, based on their relative location in respect to the closest virtual outcrop (Fig. 4). Subarea A was treated differently, since the kinematics of the Pirola fault highly affect the fracture orientations especially in the northern sector. Therefore, distinctive fracture sets from other sites and virtual outcrops are observable, as a consequence of the presence of a major shear zone (Romano et al. 2020). Subareas B and D show a distribution of fracture sets similar between data from virtual outcrops and measuring sites. In the case of VO3, an additional cluster was recognized from the virtual outcrop only (i.e., cluster 5 in Fig. 4, subarea D).
In the subarea C, two clusters of joints were measured in the field only, while a very limited number of fractures with this angle were observed in virtual outcrops and were not recognized as a cluster (i.e., clusters 3 and 4 in Fig. 4, subarea C).

Length and aperture measurements
Fracture dimensions (e.g., fracture length in Fig. 5a) of fieldbased and UAV-based measurements were compared. Both the datasets show a power-law distribution, with a possible upper bound for fracture length distribution which can be recognized at around 100 m. A systematic underestimation in the field-based data was presumed, due to the limited exposure, discontinuity of outcropping strands and nonvisibility of the fractures. The underestimation of fracture length by field measurements resulted in extremely low permeabilities (i.e., ca. 10 −24 m 2 ), if assumed as valid, presumably caused by the extremely low fracture connectivity inside the rock volume. Data comparison thus suggests an upscaling factor of the field observations by one order of magnitude, which is generally proposed for dimensions upscaling when measures are performed at outcrop scale (Chesnaux et al. 2009; Voeckler  Fig. 4 a Map showing the four regions selected for comparison between field and virtual outcrops (A-D, indicated with Thiessen polygons of different colors) and b their respective stereonet plots (poles to planes). All data collected from the field are reported as circles in the stereoplots, while all data collected from virtual outcrops are reported as crosses. The reader is referred to Table S2 of ESM1 and Figs. S1-S4 of ESM2 for details of the different recognized clusters and Allen 2012). The data used as input for the model, if coming from field data alone, were consequently increased by one order of magnitude.
Additionally, as observed from UAV-based fractures, aspect shape ratio results (Tables S1 and S2 of ESM1) are not reliable based on field observations, due to the limited exposure of fractures in outcrops. Therefore, for only field-based data, the aspect ratio input for DFN model generation was assumed to be the average value obtained from UAV measurements (Table 1; Bonnet et al. 2001).
Even if some studies (e.g. Bisdom et al. 2017b) suggest the use of outcrop network geometry for deterministic flow modeling, in the case of 3D DFN stochastic models (as in this proposed case), the need for fracture aperture, even for a qualitative flow evaluation result, is necessary, since this value highly affects permeability outputs (Seers and Hodgetts 2014). In this sense, this is a major limitation of UAV measurement for DFN modeling of water flow. Following the aforementioned considerations on the correct upscaling of fracture dimensions, and assuming the data clustering described previously, the aperture/length scaling of the fractures was compared (Fig. 5b) by using field-based aperture values and the associated UAV-based length, if available. If only field-based values were available, the upscaled value for fracture length was used. The dataset from this study fits well with other published data (a dataset spanning a length range from 10 −2 to 10 4 m) and confirms that most of the fracture clusters show a roughly square root dependence between dimension and aperture, in accordance with other studies (Klimczak et al. 2010;Seers and Hodgetts 2014).
In this study, when UAV measures only were available (and therefore aperture data input would be missing), a root square ratio of fracture length was applied, starting from a minimum value of 1 mm, presenting the 5 th percentile of the aperture measured on site.

Integrated model construction and output
Following the abovementioned considerations, the integrated model was performed as follows.
All the Thiessen polygons in the fault area (i.e., subarea A in Fig. 4) were processed using field data only, since UAVbased measurements did not show accordance with all the measuring sites in the fault zone. In the other areas data from field and virtual outcrops were integrated.
The main results from the DFN modeling are shown in Fig. 6: permeability has been converted to hydraulic conductivity, and the directions of maximum horizontal conductivity and its magnitude are reported for each of the Thiessen polygons (for a full report of the results see Table S3 of ESM1). The higher values of permeability were obtained in the SE sector of the study area, reaching a maximum of ca 1.8 × 10 −7 m 2 , equivalent to a hydraulic conductivity of 1.35 m/s with an almost isotropic ellipsoid for the polygon FO18 (Fig.  6) and a strongly oblate one for polygon FO20, with a preferential conductivity NE-SW oriented. Such a high value of conductivity results from the presence of open fractures at surface, with apertures up to centimeters. The boundary condition for a laminar flow in fractures could not be honored locally, so this value was considered as an upper bound for expected conductivities. Nonetheless, the two aforementioned polygons can potentially contribute as source areas for the springs V07, V09 and V10 (Fig. 6): the only ones showing high and constant values of water flow (up to 300 L/s for V07) for the whole summer season.
The lowest values were calculated for the area located north of the Pirola fault, with an equivalent permeability just a few orders of magnitude greater than the primary permeability of Fig. 5 a Normalized cumulative distribution for UAV-based and fieldbased fracture measurements. b Aperture/length scaling of measured fractures in the study area compared with a dataset of opening-mode fractures (after Klimczak et al. 2010); fracture lengths scale to fracture apertures as power-law functions with exponents n ∼ 0.5 the bedrock, (i.e., 10 −12 m 2 equivalent to a conductivity of 10 −6 m/s). This strong compartmentation of the aquifer is ascribable to the Pirola fault, that also results in a subparallel, ca. E-W-striking, fracturing of the bedrock in its surroundings, and, in turn, results in strongly oblate and subparallel ellipsoids for the nearby polygons (e.g., FO7 in Fig. 6).
The area immediately south of the lake and the SW area present intermediate permeability values (with an average value of 2.6 × 10 −8 m 2 ) and a preferred direction of hydraulic conductivity (1.6 × 10 −1 m/s on average) oriented ca. NE-SW. Nonetheless, two polygons (e.g., FO21 in Fig. 6), show considerably low values of conductivity, with an ellipsoid strongly oblate with K1 and K2 lying on an almost horizontal plane. These two sectors thus potentially constitute a groundwater divide, almost corresponding to the hydrologic one, separating the NE area, from the SW one.

Water chemistry results and conceptual model of water flow
Water chemistry results: selection of the best geochemical markers Water chemistry in the study area revealed, as expected, a general low mineralization: the maximum value of EC is 98 μS/cm (details of chemical data are presented in Table S4 of ESM1). Observing the Piper plot (Fig. 7a), the chemical facies of all the sampled water is Ca-Mg-HCO 3 type, in accordance with other studies on mafic and ultramafic bedrock (Harmon et al. 2009;McClain and Maher 2016;Binda et al. 2018). The main observed differences among springs are present in the cations plot, especially for Ca and Mg. This is an index of the geochemical differences between the lithologies separated by the Pirola fault, with a higher Mg load caused by serpentinite group minerals. Comparing their relative abundance, in fact, these ions present a clear inverse correlation in the relative abundance.
Nonetheless, principal component analysis better elucidates the marking of both major ions and minor metals analyzed. The loadings of variables representative of gabbros and gneiss minerals, i.e., Ca, Na, K (Zhu et al. 2003) are well separated from variables presenting a higher load for elements generally derived from the leaching of serpentinite minerals, i.e., Cr, Ni and Fe (Heikkinen and Räisänen 2008), observing the PCA plot of components 1 and 2, which summed explain the 42% of the total variance (Fig. 7c). The former present negative values of both components 1 and 2, while the latter present positive values, perfectly remarking the geochemical difference in element dissolution of springs in the study area. Therefore, while a partial clustering of different springs is  Fig. 7 a Piper plot for all water sampling points, generated using average values of major ions in all sampling campaigns. b Detail of the Ca-Mg relative abundance scatterplot. c Loading plot of principal components 1 and 2, generated with all the data collected in all the sampling campaigns observable in the Ca-Mg plot in Fig. 7b, PCA better highlights the fingerprint of all elements typical of serpentinite rocks and synthetizes the main differences in geochemistry already in principal components 1 and 2. In this plot (Fig. 7c), Mg does not correlate with Ca and other trace elements; instead it shows a correlation with EC and carbonates, because this element shows a general high concentration in all the water samples. This case study confirms the potential of PCA for geochemical analysis of water samples (Binda et al. 2020b;Gambillara et al. 2013;Giussani et al. 2016;Rosen et al. 2018).

Comparison of geochemical markers and modeled conductivity
The average scores for the principal components 1 and 2 (PC1 and PC2) are interpolated in Fig. 8 using ordinary co-kriging as an interpolator algorithm, to elucidate the differences in geochemical markers in the analyzed springs, and to compare them with the modeled flow (Gambillara et al. 2013;Hilberg 2016). The overall permeability and the direction of K1, obtained from tensor ellipsoids and projected onto the horizontal plane, are interpolated using a spline algorithm and following the procedure after Gumiaux et al. (2003).
No direct measurements are available for the piezometric surface, but groundwater flow can be presumed toward the known springs and following the main conductivity directions (Welch and Allen 2014). Starting from this assumption, and through the comparison with K1 directions, two main aquifers can be identified (Fig. 8). The Pirola aquifer is draining groundwater toward the NW, feeding the Pirola lake and an aligned series of springs along the Pirola fault. The Ventina aquifer drains toward the W and SW and is characterized by increasing conductivity values to the south; waters percolating in the southern sector of the Pirola hydrological basin are partly drained in the Ventina aquifer (see white arrows in Fig. 8). The proposed flow path is also consistent with the chemical similarities between the springs on the right side of Ventina Valley (V09 and V10) and the P08 spring.
Chemical data, focusing on the information given by scores of principal components 1 and 2 (Fig. 8), divide the study area into two main clusters of springs: (1) the spring at the Pirola Lake inlet (P04, P05 and P06) and P08, showing with low values of scores, and (2) the springs fed by the Ventina aquifer, along the valley slopes. The two aquifers are partially separated by conductivity contrasts, to the north, and by diverging conductivity directions, to the south.
In the western area, the progressive increase of values for principal components 1 and 2, from P15 to P12 and P13 ( (scaled black bars). The trace of a groundwater divide, separating the fractured aquifer into two subvolumes, is indicated (thick blue line), along with suggested water flow path (white arrows) between the waters from the fractured aquifer with the shallower phreatic aquifer, mainly hosting waters flowing through talus deposits. The mixing of the waters from the two aquifers happens in correspondence with a lateral moraine, ascribable to the Last Glacial Maximum (Fig. 1), which works as an aquitard in the shallow phreatic aquifer, causing water stagnation and promoting infiltration where a concurrent lowering of conductivity values in the fractured aquifer causes a concurrent raise of water table from that source (Hilberg and Riepler 2016).

Approach applicability and suggested workflow
The application of this integrated approach opens the possibility for a fast and cost-effective preliminary flow estimation in remote areas, where typical direct measurement and classical measuring systems are not easily applicable. In this context, the proposed approach for a complete and fast evaluation of water flow of fractured aquifers in remote settings would include a preliminary UAV flight (to recognize the main clusters of fractures) followed by a limited measurement campaign of recognition of the same clusters observed by UAV and their aperture measurements, avoiding long sampling campaigns and reducing the necessary number of measuring points (Bisdom et al. 2017a;Akara et al. 2020). Comparison of data for a conceptualization of the basin could be easily provided by geochemical markers, as well as isotopic ones (Vallet et al. 2015;Hilberg 2016;Achtziger-Zupančič et al. 2017a), requiring only periodical sampling of the springs. Nonetheless, the main drawback of this approach is the difficulty in reaching a quantitative understanding of groundwater flow through numerical detailed modeling, requiring validation through direct flow data (Massaro et al. 2018;Romano et al. 2020).
However, the proposed approach results are easily applicable in high-mountain settings and other areas presenting complex tectonic settings, where the local stress field is not easily measurable. As opposed to other studies in the literature (Agosta et al. 2010;Vollgger and Cruden 2016), this study elucidates the potential for a fast and objective preliminary observation of the main cluster of fractures and their effects on water flow, without assumptions derived from the stress field, permitting the research project to reduce time and costs for hydrogeological analysis.

Conclusions
This work aims to compare the directions of the main hydraulic conductivity components in a fractured aquifer, calculated from integrated field and UAV-based fracture measurements, with a geostatistical analysis of geochemical markers analyzed in spring water chemistry.
A comparison of UAV-and field-data collection methods detected a general accordance in fracture orientation and a more reliable assessment of fracture dimensions and aspect from UAV data. Nonetheless, an important feature in hydrogeological studies-fracture aperture-is not measurable using UAV-based data alone. The adoption of fieldbased aperture measurements for corresponding clusters of virtual outcrops, and the space discretization by means of Thiessen polygons, was revealed as a valid approach for DFN modeling. Possible groundwater flow paths have been recognized based on the distribution of known springs, and the extents of two partially separated aquifers have been mapped.
Finally, whilst still requiring on-field analysis, UAV-based measurement of fractures for water flow modeling can be a helpful, fast and cost-effective tool for preliminary exploration of fracture clusters, in order to combine orientation and dimension data of fractures with field-collected fracture-opening measurements. This integrated approach reveals a great potential of application, especially in remote settings, where detailed hydrogeological analyses are hardly achievable.
Supplementary Information The online version contains supplementary material available at https://doi.org/10.1007/s10040-021-02304-4. measurements and fractured aquifer characterization'). The fractures are regarded as a series of penny-shaped cracks, with their geometry varying in function of fracture orientation, fracture trace length, and fracture aperture.
The calculation of the permeability tensor that the Fracture Modelling module performs is based on the principles described by Oda (1985) which is itself based on the equation by Snow (1969), describing the laminar flow of fluids along smooth fracture systems: where K i is the conductivity along the ith fracture system, of frequency f i (m −1 ) and aperture e, and ν is the water kinematic viscosity. As the calculation is based on the laminar flow model between parallel plates, which includes the cube of the fracture aperture, it is very sensitive to input aperture size (e.g. Witherspoon et al. 1980). The software sums up the contributions from all the fractures included in the volume (using Eqs. 13 and 15 in Oda 1985) and gives a second order tensor, oriented in a preset reference system. Calculating the eigenvalues and the eigenvectors of the thus obtained matrix, the three main permeabilities, k1 (or k max ), k2, k3 (or k min ), and their orientations in space are obtained. With these results, it is then possible to build the permeability ellipsoid having k1, k2, k3 as semiaxes, expressing the anisotropy of the fractured rock mass, considered as a continuum medium. To visualize orientation and degree of anisotropy of permeability tensors, ellipses and ellipsoids can be used in two and three dimensions, respectively. The resulting permeability tensor is defined in the model coordinates. The eigenvalues of each ellipsoid were normalized on its k1 and used then to represent the local anisotropy of permeability. For more details about solving the permeability matrix, the reader is referred to the work by Lang et al. (2014).
Finally, if needed for discussion, permeability was converted into hydraulic conductivity by assuming a water density of 1,000 kg/m 3 and a temperature of 10°C, thus obtaining a conversion factor from permeability to conductivity of 7.5 × 10 6 (Feseker 2007).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.