## 2 Materials and Methods  
The groundwater flow was simulated for an unconfined aquifer in a steady-state condition, focusing on the saturated aquifer zone according to the Darcy law [19], through which the filtration velocity, \(v\), can be expressed by Equation (1):  
\[v=-K\times\frac{\partial h}{\partial s}=-K\times i \tag{1}\]  
where \(K\) is the hydraulic conductivity of the medium, and \(i\) is the hydraulic gradient, expressed as the ratio between the hydraulic head difference, \(\partial h\), and the spatial coordinate of the filtration direction, \(\partial s\). As the Darcian velocity is generally not available but can be estimated, the aim of mathematical-physical models is properly to assess the hydraulic head for different scenarios and stress conditions. Considering an incompressible fluid, combining the mass conservation law with the Darcy equation, a partial derivative equation was obtained, describing the flow yield of groundwater in a steady state condition [20; 21]:  
\[\frac{\partial}{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+ \frac{\partial}{\partial y}\left(K_{yy}\ \frac{\partial h}{\partial y}\right)+\frac{\partial}{ \partial z}\left(K_{zz}\ \frac{\partial h}{\partial z}\right)=W \tag{2}\]  
By using Equation (2), the water extracted for the volume unit from the porous medium (\(W\)) based on the coordinate of the filtration direction, \(\partial x\), \(\partial y\), \(\partial z\), can be estimated, considering the principal component of the hydraulic conductivity tensors \(K_{\text{xx}}\), \(K_{\text{yy}}\) and \(K_{\text{zz}}\), aligned with the x-, y-, and z-axes [22]. The hydraulic approach developed in order to analyze and model groundwater was described for a generic aquifer totally drained by a unique spring outlet, and was then applied to the Caposele spring catchment. The general model can be considered a benchmark that is then applicable to aquifers with the same hydraulic and hydrogeological properties as for the study area. The MODFLOW numerical code was used, simulating groundwater flow in the saturated zone of the aquifer, based on the conceptual model designed for the study area [4]; then, the code was implemented through the graphic user interface PROCESSING MODFLOW 5.3, which simplified the calculation code. The MODPATH package [23] allows one to calculate the particle paths and travel times of the particle trajectory in order to define the flow net, using a semi-analytical algorithm to draw the water particle. Figure 1 highlights a plan view (horizontal) of the groundwater model of a generic aquifer (obtained by 3D flow simulation) drained by a unique spring outlet and characterized by a regular shape. The groundwater modelling was carried out simplifying the groundwater framework [24], where physical and hydraulic parameters were mathematically integrated on the vertical plane, as equipotential lines, [25; 26; 27; 28; 29]. Figure 1 shows the flow net obtained for an unconfined aquifer, square shape with a 7500 m length of side, grid characterized by 22.500 cells, and with a dimension of \(50\times 50\) m; vertically, the model discretization was carried out for a maximum model depth of 2000 m considering two equally sized layers. In accordance with the MODFLOW simulation, a regular flow net was drawn, in which the flow lines converge towards a unique spring outlet. Assigning specific boundary conditions to the model, the hydraulic head of the spring was considered constant during the flow simulation (_Dirichlet type condition_), and the recharge amount was modulated on the spring discharge (_Neumann type condition_). In this model, the hydraulic head progressively decreases towards spring; conversely, the slope of the water table progressively increases towards spring. Even if the flow net shape depends on the geometry of the aquifer, the recharge amount and the hydraulic conductivity of the medium rule the maximum hydraulic head in the aquifer, following which the hydraulic head value drops between two consecutive equipotential lines.