A Double Scale Methodology to Investigate Flow in Karst Fractured Media via Numerical Analysis: The Cassino Plain Case Study (Central Apennine, Italy)

A methodology to evaluate the hydraulic conductivity of the karst media at a regional scale has been proposed, combining pumping tests and the hydrostructural approach, evaluating the hydraulic conductivity of fractured rocks at the block scale. Obtaining hydraulic conductivity values, calibrated at a regional scale, a numerical flow model of the Cassino area has been developed, to validate the methodology and investigate the ambiguity, related to a nonunique hydrogeological conceptual model. The Cassino plain is an intermontane basin with outstanding groundwater resources. The plain is surrounded by karst hydrostructures that feed the Gari Springs and Peccia Springs. Since the 1970s, the study area was the object of detailed investigations with an exceptional density of water-wells and piezometers, representing one of the most important karst study-sites in central-southern Italy. Application of the proposed methodology investigates the hydraulic conductivity tensor at local and regional scales, reawakening geological and hydrogeological issues of a crucial area and tackling the limits of the continuum modelling in karst media.


Introduction
Carbonate rocks have a heterogeneous and anisotropic hydraulic conductivity difficult to characterize, with a variation in space and in geological time [1]. The hydraulic parameters depend on the hierarchal rock structure with a facilitated circulation in the downgradient direction [2]. With a negligible permeability matrix, interruptions in the continuity of the material (i.e., faults, fractures, dissolution conduits, and coarse-grained channel fills) constitute the preferential path for the groundwater discharge and contaminant transport. For their economic importance [3], karst media have been fully investigated, including the definition of water balance [e.g., [4]], use of tracers [5], correlation between the rainfall rate and the spring hydrographs [6], hydrochemical samplings [7], or mixing hydrogeological and geophysical data [8]. Usually, followed approaches take into account mean hydraulic parameters, valid at regional scales. Modern methods investigated the fractured and karst media by a discrete approach: Atkinson [9] described the groundwater flow in a Carboniferous karst aquifer combining turbulent conduit flow and Darcian flow in fine fractures; Andersson and Dverstorp [10] predicted the groundwater flows through a network of discrete fractures statistically generated, via a Monte-Carlo simulation; Berkowitz [11] analysed open questions of flow and transport in fractured geological media; Maramathas and Boudouvis [12] described a deterministic mathematical method by the characterization of the fractal dimension of the network; Pardo-Igúzquiza et al. [13] adopted a method to generate conduit geometries via a statistical approach, based on speleological surveys; Borghi et al. [14] proposed a stochastic method to generate a karst network; Bauer et al. [15] studied the karstification process of a single conduit via a hybrid continuum-discrete approach; Liedl et al. [16] used a continuum pipe-model describing the evolution of a karst media; Langevin [17]  (1)  [4,[23][24][25]42]. Key to the legend: (1) karst springs; (2) karst streambed seepages; (3) groundwater flow direction through the karst hydrostructure.
flow model describing the fracture zone; Masciopinto et al. [18] focused on the fractures and a laminar/nonlaminar flow. However, improvements of the porous equivalent approach have been adopted in fractured or karst media, via numerical simulation of the flow [19,20] or transport model [21]. Karst media offer several solutions that can be applied at the time, involving a discrete approach, a porous equivalent, or a combination of both [22]. In this work, we investigated the fractured karst rocks of the Cassino plain, located in Central Italy ( Figure 1) and fully studied since the 1970s, due to its outstanding groundwater resources. Celico [23][24][25] and CASMEZ [26] focused on the Mesozoic and Cenozoic karst layers with a high density of pumping tests and other indirect investigations. Nevertheless, geological and hydrogeological setting of the Cassino plain is debated with two different conceptual models ( Figure 2). Area is appropriate for a continuum approach in the study of the fractured karst rocks at the regional scale. In the area, investigated karst layers were not affected by the huge karst processes related to the Mediterranean salinity crisis event [27], due to the recent uplift of the Apennine orogen [28]. The subsequent normal faulting event deposited thick continental sequences in the low-stand sectors, with a function of permeability boundary. These sequences affected the groundwater circulation, hindering the development of a karst network in the spring areas. The combination between tectonic activity and sedimentation rate varied the base level of fractured carbonate aquifers [29]. As a consequence, basal springs are outstandingly high and steady, without impulsive response to seasonal recharge. In accordance with other hydrogeological Venafro Mts.  [26]. No karst caves have been reported in the original documents; (c) redrawn hydrogeological model from Boni and Bono [42]. The authors [4,42] schematize the carbonate bedrock via a horst and graben setting and a flux from Gari Springs (G1, G2, and G3) to Peccia Springs. (d) modified hydrogeological section from Celico [23]. Celico [23][24][25] supposed the interruption of the carbonate bedrock and no underground connection between Gari Spring and Peccia Springs.  Geomechanical analysis including the evaluation of orientation, spacing, and aperture (a) allows constructing discrete-fractured cubes (b) representing the average fracturing condition. Numerical models at 1 m scale (c) have been calibrated via pumping tests (d) estimating the hydraulic conductivity tensor at regional scale (e). The latter has been tested via numerical simulation at regional scale and compared with conceptual models from literature (f). studies in Central Apennine [30], groundwater flows through fracture networks, instead of karst conduits. Thus, fractured karst masses of the studied area cannot be compared with the other Mediterranean karst basins or older Paleozoic carbonate rocks.

N CWW
To characterize the groundwater flow through the carbonate layers of the studied area and verify the importance of the fracture networks over the karst conduits, a methodology has been proposed to estimate the hydraulic conductivity tensor ( Figure 3). This methodology starts with geomechanical surveys to define average fracturing setting parameter, further implemented by a FEM (Finite Element Model) analysis at a 1 m scale. This first step is in accordance with the hydrostructural approach [31][32][33]. Estimated dataset has been calibrated with pumping tests to derive a reliable hydraulic conductivity tensor and successively validated with a Darcian numerical analysis at a regional scale. Methodology, validated for the Cassino plain, can be applied across the Central Apennine, thanks to the similar hydrogeological setting. Based on the new collected data and the review of the literature, this study aims to (i) define a hydraulic conductivity tensor of the investigated fractured karst media, representative at regional scale; (ii) verify the adoption of a downscaling procedure, able to shift from 1 m scale model to a regional one; (iii) offer an in-depth analysis related to the highly debated conceptual models of the study area.

Geological and Hydrogeological
Setting. The Cassino plain is an intermontane basin, located downstream of the Latina Valley in the Central Apennine ( Figure 1). Mesozoic dolostones constitute the base of the stratigraphic sequence, followed by Mesozoic and Cenozoic limestones. The entire carbonate sequence has a thickness of 2000 m [34]. Upwards in the sequence, Tortonian sandstone and clay levels represent the foredeep deposits of the Apennine chain, with a thickness of 800 m [35]. Quaternary continental deposits follow [36]. The spatial distribution of the stratigraphic sequence is strongly controlled by tectonics. Since the Messinian [28,37], the Mesozoic and Cenozoic carbonate deposits experienced contractional tectonics and were structured in a thrust and fold geometry. The latter was completely masked by strikeslip tectonics along the Val Roveto-Atina-Caserta line [38,39] and a subsequent normal faulting event. The last tectonic event dissected the Mesozoic and Cenozoic carbonate sequences by normal faults with NE-SW and NW-SE trends, creating a horst and graben geometry. Indeed, the flat alluvial morphology of the plain is interrupted by two main monoclines (Trocchio Mt. and Porchio Mt.) made of Mesozoic and Cenozoic limestones and bounded along three sides by normal faults (Figure 2(a)). Numerous minor carbonate structures outcrop from the alluvial plain. Unfortunately, in literature a unique structural setting has not been agreed: Zalaffi [40] assumed Trocchio and Porchio Mts. as horsts but, two years later, Accordi interpreted the latter as klippe [41]. The problem related to the structural setting was accentuated Geofluids 5 in case of hydrogeological modelling [4,[23][24][25]42]. All the conceptual models are in accordance with defining the Gari Spring (elevation of 41-30 m a.s.l., discharge of 18 m 3 /s) and Peccia Spring (elevation of 29-25 m a.s.l., discharge of 5 m 3 /s) as related to the karst Mesozoic-Cenozoic aquifers. From a hydrogeological point of view, the sequence can be schematized in four units. The dolostone unit (DU) represents the basal aquifers (Figure 2(a)). Though it has several fractures and karst elements, the hydrogeological characteristics are similar to terrains with interstitial hydraulic conductivity and high storage capacity [4]. The Cinquina water-well (CWW in Figure 2(a)) suggests a hydraulic conductivity of 5 ⋅ 10 − 5 [m/s]. Conversely, a continuous network of fractures affects the younger limestone unit (CU), contributing to a higher conductivity value. The limestone unit was investigated via 36 boreholes furnished with pumping tests [26], during the construction of the western Campania aqueduct (WCAWW in Figure 2 In the end, the quaternary unit (QU) has a variable importance related to their extension ( Figure 2(a)). At a regional scale, these deposits act as an aquitard of the limestone and dolostone sequence. From the reconstructed hydrogeological sequence the importance of the limestone unit (CU) and its spatial characterization is evident. Boni and Bono [42] linked the Gari Springs to the limestone sequences. They individuated in the Gari Spring area (Figure 2(a)) several carbonate horsts, delimitated by the quaternary unit. This setting gives origin to localized springs with different elevation: the Cassino urban area springs (G1 in Figure 2

Hydraulic Conductivity Evaluation.
The study of the groundwater flow through the karst media requires the definition of the hydrodynamic parameters of the carbonate layers, their spatial variability, and the scale of the analysis. At regional scale, the hydraulic conductivity is one of the most important parameters and its definition allows characterizing the groundwater flows as well as the contaminant transport processes. Because groundwater flows along preferential ways (e.g., fractures and conduits), the magnitude of the hydraulic conductivity should be represented along the three main directions, for full characterization. The Cassino plain (Figure 1) represents an ideal case for the study of the fracturing setting; the carbonate layers outcrop along the whole plain and the area was fully investigated from a geological and hydrogeological point of view. To estimate the hydraulic conductivity in fractured karst media, we calibrated a methodology combining the hydrostructural approach and the pumping test data ( Figure 3). The hydrostructural approach [31][32][33] requires geomechanical surveys to recognize and characterize the joint sets (i.e., a group of subparallel joints) (Figure 3(a)). The approach starts with geomechanical measures, whereas the karst layers outcrop. A total number of 40 geomechanical stations have been performed. Every geomechanical site is located where the bedrock is well exposed, without detritus or vegetation. If possible, measurement sites have been chosen in correspondence with excavations or quarries. Indeed, two exposed rock-walls with different orientations better describe the three-dimensionality of the rock masses. In addition, geomechanical stations were also located nearby roads or other linear features. According to the principles of the geomechanical standards [45], the joint sets have been identified, taking note of their spacing and opening. For our purposes, the joints were assumed smooth and plain; indeed the roughness of the joint wall along the topographic surface is not representative of the rock-mass condition in depth. Because the fracturing condition has a strong spatial variability difficult to characterize, in terms of density and characteristic of the joints, the definition of a mean fractured setting of the karst media was necessary. After the definition of the recognized joint sets, different values of spacing and opening have been considered, constructing a discrete fracture network (Figure 3(b)). Thus, the different fracturing cases have been considered. Fractures have been bounded in cubes with a side length of 1 m, oriented with the following annotation: positive = east, positive = north, and positive = up. For every fracturing setting, numerical simulations by a FEM code have been performed to calculate the flow through the faces of the cubes (Figure 3(c)). A FEM gives the opportunity to easily modify the opening of the fractures, not varying the geometry of the model. Due to the high heterogeneities of the limestone layers, the spacing dataset has been separated in three intervals as well as the opening dataset in four intervals. The combination of the two parameters (spacing and opening) gave back a final dataset of twelve different fracturing cases. COMSOL5 software and its Fracture Flow Interface in the Porous Media and Subsurface Package were used. Inside the single fractures, the following equations have been assumed: where ∇ is gradient along the fracture extent; is average fracture opening; is flow velocity; is liquid density; is  outflow; is fracture permeability; is fluid viscosity; is pressure.
According to the cubic law [46], where ℎ represents the hydraulic head variation and a form factor; the permeability along a single fracture ( ) can be written in the following form [47]: Though several studies outlining the limit of the cubic law [48,49], it is a point of reference for the rock mechanics and it is usually used for flow in fractured systems [50][51][52].
During the construction of the numerical model, fractures were represented as planes subdivided via a freetriangular mesh. The number of the triangles varies in function of the considered spacing set as well as the number of fractures inside the single cube; usually it is around 60,000. The dimension of the triangles strongly diminishes near the intersection between two planes. The flow through the rocky matrix was neglected, assuming that water cannot flow outside the fracture discrete network (Neumann condition type). As supplementary boundary condition, an average hydraulic gradient of 0.5%, typical for the limestones of the Central Apennine [4,25] has been applied (Dirichlet condition type). The gradient was applied firstly along the -axis. Simulations were repeated applying the gradient along the -axis and the -axis. Measuring the simulated discharge, a porous equivalent hydraulic conductivity has been calculated for every cube and along the three directions. Therefore, the hydraulic conductivity dataset was computed, varying the fracturing condition. Dataset has been calibrated with pumping tests (Figure 3(d)), in order to individuate a permeability tensor (Figure 3(e)) valid at regional scale. Numerical analysis of the flow model at regional scale (Figure 3(f)) allows validating the proposed methodology for the study of karst media.

Geometry and Boundary Condition at Regional Scale.
The calibrated hydraulic conductivity tensor has been validated via numerical analyses at regional scale. A geological section for numerical purposes has been constructed ( Figure 4) along the A-B section trace of Figure 2(a), assuming the horst and graben conceptual model [4,42], with some adjustments due to new surveys and the reinterpretations of borehole stratigraphies [25,44]. Because the Cassino plain has a nonunique geological and hydrogeological setting, the most recent conceptual model [4] has been preferred. The model has a strong horizontal variability due to the presence of normal faults; therefore a FEM analysis has been preferred, using COMSOL5 software and its Porous and Subsurface Package. According to the literature review, 40 domains have been classified in four hydrogeological units: the dolostone unit (DU), the limestone unit (CU), the terrigenous unit (TU), and the Quaternary unit (QU). At regional scale, every unit is considered homogeneous and porous and the Darcy Law valid [53]. The entire geometry was discretized via 39.457 elements, in a free-triangular mesh, with refinements in correspondence with the domain vertices. Because the numerical model just verifies the final part of the hydrostructures, a steady-state hydraulic head, from the inner Cairo Mt. to Gari Spring and from Camino Mt. (Venafro Mts.) to Peccia Spring, was applied (Dirichlet condition type). The hydraulic gradient of 0.5% has been chosen in accordance with Boni et al. [4] and Celico [25]. Along the vertical borders and at the bottom, no flow conditions (Neumann condition type) have been applied. To simulate the unsaturated aquifer thickness, a Heaviside function of the pressure (value from 0 to 1) has been considered along the topographic surface generating a hydraulic head dependent boundary (Dirichlet type-condition), according to Chui and Freyberg [54]. The latter condition cancels possible undesired fluxes towards the unsaturated zone, very common during the numerical analysis. Boundary conditions are represented in Figure 4.

Hydraulic Conductivity Dataset from the Hydrostructural
Approach. The geomechanical analysis of the rock masses pointed out a wide dataset of fractures, collected in correspondence with the dolostone and limestone outcrops. The dolostones are located at the northeast boundary of the Cassino plain between Cervaro and San Vittore (Figure 2(a)). In this area, the dolostone layers are tectonically disturbed and weathered. During the surveys, a systematic geomechanical analysis of the fractures was not possible. The fracture density is very high and the rock-mass has a low compactness, assuming characteristic of a loose, cemented sand. Surveys pointed out fine silt in the main fractures. The hydraulic conductivity measure of the dolostone ( = 10 − 4 m/s) coming from the Cinquina water-well (Figure 2(a)) is in accordance with the geomechanical observation. Zalaffi [55] described the same hydrostructural setting of the Matese Mts. dolostones ( Figure 1).
Conversely, the hydrostructural approach has been fully applied to the limestone layers, between Gari and Peccia Springs. The results of 36 geomechanical stations pointed out a well-organized fracturing setting with several joint sets that take origin from the stratigraphic and tectonic features of the limestone rock-mass. The dip-slip measures of the bedding vary from N.90 ∘ to N.130 ∘ across the studied area. Data are in accordance with the schematization of Boni et al. [4] and Celico [25]. The dip measures are averagely around 30 ∘ except at the Montecassino Hill, where a recumbent anticline fold structure induces a local variation of the dip between 10 ∘ and 80 ∘ , in accordance with the geological data of Saroli et al. [44]. The spacing (i.e., thickness of every strata) is averagely around 15-40 cm but locally can reach values of 1 m. The opening of these joints is very low, around 0.1 mm, but it is locally higher when the fractures are affected by flexural slip. The bedding can reach an opening of 1 mm. The tectonic joints instead have not a homogeneous dipslip, and a schematization in different joint sets is required. Collected data shows a diffuse high-angle dip, from 70 ∘ to 90 ∘ . According to the geomechanical standards [45], all the fractures, with high-angle dip-slip, have been gathered in the same group isolating four tectonic joint sets: the fractures with dip angle in the interval between N.345-N. 15  . The spacing and the opening of these tectonic joint sets are variable and depend on the local heterogeneities that affect the rock-mass. After the individuation of the joint sets, the spacing and opening dataset have been divided into classes, to reproduce different degrees of the fracturing condition: in total three spacing classes and four opening classes. A summary of the collected data is in Table 1.
Assuming the simultaneous presence of all the joint sets, cubes with a discrete fracture network have been constructed. The orientation of the fracture has been considered constant, assuming an average value. Then, the three spacing sets have been crossed with the four opening-sets, for a total of twelve numerical analyses. According to the boundary condition, simulated fluxes occur only through the planes that represent the discrete fractures. The numerical simulation pointed out a homogeneous dissipation of the imposed hydraulic head across the opposite face of cubes along the -axis ( Figure 5(a)), -axis ( Figure 5(b)), and -axis ( Figure 5(c)).
Black arrows ( Figure 5) outline how the simulated flux is parallel to the fracture orientation, respecting the boundary condition. The ratio between the discharge along the selected direction and the imposed hydraulic head is proportional to the hydraulic conductivity value, if an equivalent porous medium is considered in accordance with the Darcy Law [49]. The simulated hydraulic conductivity range is very wide, from 10 − 12 to 10 − 3 [m/s] (Figure 6). According to the model, the hydraulic conductivity strongly varies in function of the chosen opening set and subordinately in function of the spacing set. Indeed, inside a single fracture each increment of the opening corresponds to a square increment of the velocity. A wide dataset of hydraulic conductivity up to nine orders of magnitude is in accordance with the huge heterogeneity of the fractured karst media.
With a wide dataset of results ( Figure 6), it would be senseless deriving a reliable hydraulic conductivity value for the limestone unit (CU), based just on the hydrostructural results. However, in every cube and are almost equal and is higher ( Figure 6) for all the spacing-opening combinations. Adding more, the ℎ/ V ratio ( + /2 ) remains approximately constant around 0.51. Thus, the results point out a value of the anisotropic rate at regional scale that pumping tests cannot provide. The 36 pumping tests show a hydraulic conductivity magnitude between 5 ⋅ − 03 and 5 ⋅ − 04 m/s, as well as the hydrostructural approach of a ℎ/ V ratio around 0.5. By the combination of the hydrostructural approach and the pumping tests, we can assume as acceptable at regional scale the following hydraulic conductivity dataset for the limestone unit (CU): CU = 5 ⋅ 10 − 4 m/s; CU = 5 ⋅ 10 − 4 m/s; CU = 10 − 3 m/s.

Limits of the Hydrostructural Approach.
The hydrostructural approach is based on the reconstruction of discrete fracture network and further numerical simulation. Geomechanical analysis reconstructs the heterogeneities of the limestone unit (CU) by punctual measures that cannot be homogeneously distributed. Indeed, vegetation or detritus mostly covers the limestone unit, especially along the piedmont band. Adding more, in absence of explorable karst conduits and caves, all the measures have been achieved just along the topographic surface. Thus, the described setting could not be representative of the fracturing condition in depth. Adding more, the cubic law [46] can be assumed valid under laminar conditions of the flux. Though in the studied area karst processes are young and not well-developed, the flux through wider fractures (karst conduits) and related turbulent conditions has been not taken into account, due to the regional character of the work and the hydrogeological setting of the Central Apennine [29,30]. However, conduits can influence the results at local scale. In this case specific laboratory tests and field-tests [48][49][50] are necessary to    investigate the relationship between fractures and conduits and the transition from laminar to turbulent flow.

Validation of the Methodology via Regional Numerical
Analysis. The estimated hydraulic conductivity dataset of the Cassino limestone unit (CU) comes from the combination between the hydrostructural approach and pumping tests results. A regional numerical simulation based on the estimated hydraulic conductivity values is a powerful instrument to validate the adopted methodology and the most recent conceptual model [4,37].  (Figure 7(a)). A similar setting is observed on the opposite side of the Cassino plain, whereas the velocity field simulates the Peccia Springs (Figure 7(a)). Nevertheless, the central part of the section shows very low velocities and a groundwater flux from Gari to Peccia Springs, through the Trocchio and Porchio Mts. (Figure 7(b)). Because the groundwater velocity field and the spring distribution perfectly matches, the next step is the validation of the conceptual model from a quantitatively point of view. In a 2D model, to verify the existence of a 4 m 3 /s through the Trocchio and Porchio Mts., we based our analysis on the simulated ratio between the downstream and the upstream Gari Spring. The latter is around 9%.
According to the Boni conceptual model [4], there is a flux of 18 m 3 /s from the Cairo Mt. to the Gari Spring and of 4 m 3 /s from Gari to Peccia Spring via the Trocchio-Porchio Mts. alignment, for a total discharge of 22 m 3 /s. Thus, the ratio between downstream and upstream discharge should be higher, around 18% at least. This ratio increases up to 30%, if we consider the minimum values of the Gari Spring discharge (13 m 3 /s). Consequently, the Boni et al. [4] conceptual model is not numerically verified, along the considered section. The numerical model takes into account favourable conditions to maximize the flux from Gari to Peccia Springs, since there are many possible solutions, compatible with hydrogeological constraints. According to Boni et al. [4], the alignment between the Gari and Peccia spring would represent a continuous groundwater flow direction similar to a conduit, laterally bounded by a thick blanket of terrigenous deposits (TU). Since no lateral groundwater exchanges occur, the 2D section model neglects possible hydraulic head loss related to the real 3D geometry. In order to reduce the hydraulic head loss from Cairo to Peccia Springs, the maximum thickness of the limestone unit (CU) has been assumed in the model. Hydraulic conductivity has been considered constant with depth up to −1500 m a.s.l., neglecting the natural permeability reduction related to the fracture closure; neither lateral variation of the permeability, related to mylonitic bands in depth, nor tectonic-tightening of the limestone unit thickness (CU) has been considered. From a hydraulic point of view, the maximum hydraulic gradient from the highest elevation of Gari Spring to the lowest of Peccia spring has been calculated. No spatial variation of the hydraulic conductivity generates a loss of the hydraulic head, disfavouring the communication between the Gari and Peccia Springs. Thus, as a precaution, a Spring. The downstream/upstream discharge ratio is around 9%, calibrated under the best assumptions, and is not sufficient to justify the Boni conceptual model [4,37].
homogeneous value of the limestone unit (CU) has been set. However, small variations of the hydraulic conductivity do not provide any significant variation on the discharge ratio. The only way to increase the upstream/downstream discharge ratio close to the theoretical value of 18% is to switch off the imposed gradient (Dirichlet condition) from Camino Mt., which have a function of backpressure to the flow from Gari to Peccia Spring (Figure 7(c)). This last condition is unacceptable and not in accordance with Celico [23][24][25] and Boni et al. [4].

Conclusions
The proposed methodology, applied to fractured karst media of Cassino plain (Figure 2), allowed defining the average hydraulic conductivity tensor at regional scale. The proposed methodology starts from the geomechanical analysis, defining the average fracturing parameters. A 1 m scale FEM analysis allowed estimating a wide hydraulic conductivity dataset, successively combined with pumping tests. Based on the literature review, application of the proposed methodology, and validation via a regional flow model, the following conclusions can be delineated: (i) 1 m scale analysis and the subsequent FEM analysis are appropriate to investigate the hydraulic conductivity anisotropy but simulated magnitude range is too wide, reflecting the heterogeneities of the fractured karst mass at local scale. On the other hand, pumping tests investigate a higher volume but assuming the media homogeneous and isotropic. (ii) Thus, a valid approach to determine the hydraulic conductivity tensor in fractured media, at regional scale, could be obtained by an accurate combination between a detailed scale analysis, based on hydrostructural scale measurements, and average values from pumping tests. (iii) The estimated hydraulic conductivity tensor is valid at regional scale, under Darcian conditions. The numerical model perfectly simulates the distribution of the karst springs along the topographic surface, according to Boni et al. [4] and Celico [23][24][25].
Nevertheless, under the best assumptions, the model is not quantitatively verified and further data are necessary to clarify the debated conceptual model of an important resource of Central Italy.
The proposed methodology represents an improvement of the continuum approach and can be applied in fractured media, especially when an average value of the hydraulic conductivity tensor is required. For fractured karst media, the Cassino plain and, in general, the Central Apennine region represents a favourable case, since the groundwater preferentially flows in the fracture networks and not in the karst conduits [30]. As future direction, methodology can be implemented with the superimposition of conduits, investigating the transition from laminar to turbulent flow. With this modification, the latter can be applied on every kind of karst media, coupling a discrete and continuum approach.