A Macroscale Hydrogeological Numerical Model of the Suio Hydrothermal System (Central Italy)

The complex behaviour of the Suio hydrothermal system (central Italy) and its potential exploitation as a renewable energy source are still unclear. To quantitatively evaluate the geothermal resource, the Suio hydrothermal system has been investigated with a hydrogeological numerical model that couples fluid flow, thermal convection, and transport of diluted species inside a hybrid continuum-discrete medium. The numerical model, calibrated and validated with available and new experimental data, unveiled the complex behaviour of the hydrothermal system. The normal tectonic displacements, the fracturing of the karst hydrostructure, and the aquitard distribution strongly influence the hydrothermal basin. In particular, a dual fluid circulation, sustained by steady-state thermal and pressure gradients, modulates the hydrothermalism at the several springs and wells. The presence of a medium to a low-temperature reservoir allows for potential exploitation of the geothermal resource.


Introduction
In recent decades, the worldwide growth of energy demand and the increase in CO 2 emissions boosted the development of new techniques for the exploitation of noncarbon sources of energy from the sun, wind, tides, and subsurface heat. In Italy, geothermal energy aroused a growing interest [1]. Indeed, the Italian geothermal potential up to economically convenient depths is considerable, with high-temperature resources (>150°C) located in the peri-Tyrrhenian sector of central Italy and in some islands of the Tyrrhenian Sea, while medium-to-low temperature resources (<150°C) are located in vast areas of the national territory [2] (inset in Figure 1(a)). Exploration and exploitation have concentrated for high and medium enthalpy fluids at shallow depth in areas of recent magmatism only [1,3,4]. However, the recent technological developments in the field have extended the potential of geothermal reservoirs to lower temperatures and greater depths [5]. Several exploration permits have been requested by private companies in Italy, indicating the significant interest of industry for this renewable resource [6]. It is therefore essential to improve the knowledge of potentially exploitable hydrothermal areas in order to increase energy production from renewable and environmentally sustainable resources.
The Suio hydrothermal basin (Figure 1(a)) shows all the characteristics of a high potential medium-to-low enthalpy geothermal system. The area shows several thermal springs with temperatures up to 50°C and gaseous emissions, located along the southeastern boundary of the Eastern Aurunci Mts, at the contact with the Roccamonfina volcanic edifice. The Suio area has been investigated by previous geological, hydrogeological, geophysical, and geochemical studies [7][8][9][10]. These works provided valuable hints about the subsurface setting of the area and allowed for the construction of a hydrogeological conceptual model of the deep and shallow groundwater flow systems. The proposed conceptual model would explain the geochemical features of the Eastern Aurunci Mts springs and the hydrothermalism of Suio basin [8]. However, none of these works validated the proposed scheme of the Suio hydrothermal basin with numerical models.In recent decades, new numerical tools of flow and transport within porous and fractured media have been developed for the investigation of hydrothermal systems. These tools allow to properly considering most of the critical features such as the lithostratigraphy, the tectonic setting, the groundwater flow, and the heat source [11][12][13][14][15][16][17][18]. The most commonly applied methods can be grouped into three main categories: (i) continuum methods, including finite difference methods (FDM), finite element methods (FEM), and boundary element methods; (ii) discrete methods, including discrete-element and discrete-fracture network methods; and (iii) hybrid continuum-discrete methods [19]. The choice between continuum and discrete methods depends mainly on the problem scale.
In this work, we developed a large-scale numerical model of the Suio hydrothermal area applying a hybrid continuum-discrete approach. The developed model has a double scope: (i) verify the conceptual model of groundwater flow and heat, taking into account the fractured nature of the system; and (ii) provide a valid tool which can be used for a quantitative evaluation of the geothermal potential of the area for future exploitation. In detail, we calculated the effect of temperature, pressure gradients, and dissolved gases on the groundwater flow inside the hydrothermal system and we verified the modelling results with available and new data. The results unveiled the complex behaviour of the Suio hydrothermal area and provided useful insights into the dynamics and exploitation of hydrothermal systems.

Data and Methods
2.1. The Suio Hydrothermal Area. The Suio hydrothermal basin is located in the southern part of Lazio region, central-southern Italy, between the Eastern Aurunci Mts and the Roccamonfina volcanic edifice (Figure 1(a)). At its southeastern edge, the carbonate complex (CC) of the Eastern Aurunci Mts is dissected and lowered by three normal faults (Figure 1(b)). The first fault has a NE-SW trend and reuses an old frontal thrust, the second has a NE-SW trend and delimits the Garigliano graben, and the third has an E-W trend and reuses an old contractional lineament [20,21]. These faults dislocate the CC complex roof to more than 1000 meters below the sea level [9] making it at contact with turbiditic (FC), sandy-conglomeratic (GC), and volcanic (VC) complexes (Figure 1(c)). For a detailed description of the geological evolution of the area, see Saroli et al. [8].
The hydrogeology of the area is driven by the vast karst hydrostructure of the Eastern Aurunci Mts with highdischarge basal springs, located at the more topographically depressed aquifer boundaries (Figure 1(b)). The other complexes have low permeability, hosting local confined, or phreatic aquifers that may feed small seasonal springs. Indeed, the Gallo well (85-1) [7] (Figure 1(a)), drilled approximately 800 m inside the Roccamonfina caldera, shows negligible fluid flow and temperatures up to 35.6°C at the bottomhole ( Figure 2). Therefore, groundwater flow entirely develops inside the CC complex, while the other complexes can be assumed as impermeable.
The conceptual model of the hydrothermal system, based on literature information, can be resumed as follows. The fluid flow is mainly oriented NW-SE into a relatively small volume (Figure 1(a)) [8]. The permeability of the CC complex is high because of tectonically induced fractures, which favours the development of karst features (e.g., solution enlargement of joints and fractures). Hydrothermal features originate from the hot metamorphic crustal rocks that delimit the bottom of the CC complex at approximately 4000-5000 meters below the sea level [7,9,22]. These rocks are hot because of the peri-Tyrrhenian volcanism [23], where the Roccamonfina volcano represents the shallowest expression. Metamorphic rocks generate heat, warming the huge aquifer hosted in the carbonate hydrostructure and feeding stationary convective loops inside the reservoir, which pull up hot fluids that subsequently mix with the cold water coming from the CC complex hydrostructure (Figure 1(c)) [8]. Indeed, the most noticeable hydrothermal effects are observable at the contact between the CC complex and the Roccamonfina volcano. Here, the mapped springs (from S2 to S19 in Figure 1(b)) and water wells (from W1 to W31 in Figure 1(b)) show temperatures and nonkarst ion concentrations that progressively increase moving towards the Roccamonfina edifice. Hydrochemical studies [7,8,10,24] show that at the NE and SW limits of

Geofluids
the Suio hydrothermal basin (S2 and S19 in Figure 1(b)), the groundwater is cold (T ≈ 16°C) with a predominance of dissolved karst ions (Ca 2+ , HCO 3 -). Conversely, moving towards the SE limit of the hydrostructure, the water temperature strongly increases up to 50°C (S9 spring in Figure 3(a)), and non-karst-origin ions (Na + , SO 4 2-, Cl -, and K + ) prevail with respect to the karst ones. Nonkarst ions cannot come from the other hydrogeological complexes bounding the reservoir because of their low permeability and minor water circulation [8]. Therefore, the recorded hydrochemical features [8] strongly suggest a continuous mixing between the cold, calcium-carbonate signature waters with hot, deep fluids. The admixed fluids have a relevant salinity due to the leaching of the karst network and possess high temperatures and non-karst-origin ions due to the interaction with the heat source hosted in the metamorphic rocks. Moreover, the hot fluids leach the karst network of the carbonate rocks and enrich in diluted gasses (mainly CO 2 and subordinately H 2 S). The rising diluted gas gradually decompresses and consequently degas, producing fumaroles, bubbling water from hot springs ( Figure 3(a)), and the supposedly artesian behaviour of water wells (Figures 2(b)-2(d)). The mean temperature of the carbonate reservoir varies between 140 and 170°C according to geothermometer estimations [8,10].

Governing Equations for Hydrothermal Systems.
To validate the conceptual model of Figure 1(c), we simulated the steady-state pore fluid diffusion by coupling the heat transfer and the fluid mass circulation with dilute species in porous media. The modelisation has been performed with the COMSOL Multiphysics ® software package [25] assuming an undeformable porous medium (i.e., uncoupled stresses and strains), laminar fluid flow regime, incompressible fluid, and ideal gases.
The heat transport equation describes the heat transport in the subsurface [26]: where ρ (kg/m 3 ) is the density of the fluid, C p (J/kg K) is the fluid heat capacity at constant pressure, ρC p eff (J/m 3 K) is the effective volumetric heat capacity at constant pressure, T (K) is the temperature, u (m/s) is the fluid velocity field, Q (W/m 3 ) is the heat source, and λ eff (W/m K) is the effective thermal conductivity of the solid-fluid system. λ eff is given by the weighted arithmetic mean of fluid and porous matrix conductivities: where λ p and λ are the thermal conductivity of the solid and fluid, respectively, and θ p (-) is the volume fraction of the solid, given by 1 − ɛ p , where ɛ p (-) is the porosity. Properly solving heat transport requires incorporating the flow field. In particular, Darcy's law can describe the fully saturated and mainly pressure-driven flow in deep geothermal strata: where the velocity field u (m/s) depends on the permeability κ (m 2 ), the fluid's dynamic viscosity μ (Pa s), and is driven by a pressure gradient p (Pa). Darcy's law is then combined with the continuity equation: where ρ (kg/m 3 ) is the fluid density and Q m (m 3 ) is a mass source term. If the simulated scenario concerns large geothermal time scales, the time dependence due to storage effects in the flow is negligible. Therefore, the first term on the left-hand side of the equation (4) vanishes because the density and the porosity can be assumed constant over time.
Fracture flow may locally dominate the flow regime in geothermal systems, such as in karst aquifer systems. Flow inside fractures are governed by the same equations governing the fluid flow inside porous media as follows: where ∇ T denotes the gradient operator restricted to the fracture's tangential plane, d f (m) is the fracture thickness, and κ f (m 2 ) and ɛ f (-) are the fracture permeability and porosity, respectively. Finally, the presence of diluted gas is modelled by simulating the transport of diluted chemical species through diffusion and convection according to the mass balance equation:    [8] has been assumed as a reference to build the numerical model. A 2D finite element model ( Figure 4) that extends 21.5 km horizontally and to a depth of approximately 6 km has been implemented. The mesh is composed of three-node triangular elements (18346 elements). The 2D geometry crosses the area characterised by the higher water circulation and hydrothermal effects, and it is almost parallel to the fluid flow [8].
The model geometry is discretized into five layers ( Figure 4), whose thickness has been estimated according to the available literature [8,9,29]: the volcanic complex (VC), the sandy-conglomeratic complex (GC), the turbidite complex (FC), the carbonate complex (CC), and the lower crust or metamorphic rock (MC). Hydraulic boundary conditions have been selected according to the hydrogeological setting of the area [8]. The dashed blue line at boundary 1 in Figure 4 simulates the mean piezometric level inside the CC complex employing a linear hydraulic head with a gradient i = 0 4% [30,31]. Boundaries 2, 3, 4, 5, and 6 in Figure 4 are impermeable (no flow orthogonal to the boundary) because the bottom and the right side of the model are bounded by impervious lithologies (MC and VC complex in Figure 1), while the left side corresponds to the watershed limit of the Eastern Aurunci hydrostructure [8]. Along the topographic surface (boundaries 7, 8, and 9 in Figure 4), a mixed boundary condition is used to split the boundary into a Dirichlet portion for the potential seepage face and a Neumann portion for the regions above the seepage face [32]. In detail, boundaries above the assumed groundwater table have a Neumann-type condition (no flow), while boundaries below the groundwater table have a Dirichlet-type condition (allowed flow).
Thermal boundary conditions consist of a heat source with constant temperature T max at the bottom left of the model (boundaries 4 and 5 in Figure 4) representing the heat flux coming from the hot metamorphic rocks close to the Roccamonfina volcano's caldera. The location and extent of the heat source have been selected according to the mean temperatures and heat flux at a depth available in the study area (inset in Figure 1(a)) [2]. The bottom right (boundary 3) and the sides of the model (boundaries 2 and 6) are assumed thermally insulated since the right side approximately corresponds to the vertical axis of the Roccamonfina volcano. Along the topographic surface, the convective heat flux is assigned to boundaries 7 and 9 with a mean heat transfer coefficient h = 5 W/m 2 K and an external mean temperature of 20°C [33], while the boundary 8 is assumed as an outlet with an orthogonal convective flow, as typical boundary conditions. Finally, boundary conditions for dilute species concentration consist of a fixed gas concentration c 0 at the boundary line 4, no mass flows orthogonally to the boundaries 2, 3, 5, and 6, and an open boundary condition is assumed for the topographic surface (boundaries 7, 8, and 9) to simulate convective outflow, with an external species concentration c 0 = 0 mol/m 3 .
The two main faults that dissect and downthrow the carbonate complex (Figure 1(c)) have been modelled as discrete fractures (the red segments a and b in Figure 4). The fracture is simplified with a set of parallel segments to study the fluid flow through it. By assuming laminar flow between the two flat parallel segments, the fluid flow obeys to the Darcy law   (5) and (6)) and the relation between fracture aperture and its corresponding permeability κ f (m 2 ) in the direction parallel to the fluid flow is derived from the wellknown cubic law [34]: where d f (m) is the fracture's aperture. The fracture permeability tensor computed with equation (9) is aligned with the local coordinate system of the fracture itself. This local coordinate system is often rotated by an angle θ with respect to the global coordinate system of the model ( Figure 4). Therefore, the permeability tensor of every fracture should be stated in a global coordinate system according to the following relation [35]: where κ f g (m 2 ) is the permeability tensor with respect to the global coordinate system and θ is the positive counterclockwise rotation angle between the global x-axis of Figure 4 and the modelled discrete fractures.

Model-Parameter Estimation.
The parameters of the different components, i.e., the fluid and solid phases, have been calibrated using the field and literature data. The simulated fluids consist of water filling the pores and CO 2 for the diluted gas, the latter being the dominant diluted species [8]. Fluid density, dynamic viscosity, specific heat capacity, and thermal conductivity vary with the temperature according to well-known experimental relationships [36,37]. For the solid matrix, each layer in Figure 4 has different properties. Mass density, specific heat, and thermal conductivity have been determined for each complex through laboratory experiments and literature data [38][39][40][41][42] (Table 1). Regarding the hydraulic parameters, i.e., permeability and porosity, the collected bibliographic data [31,43] outline that fluid diffusion develops essentially inside the carbonate complex (CC), which is permeable by fractures and karst processes, while the other lithologies can be assumed as impervious [8,43,44]. Therefore, the VC, GC, FC, and MC complexes present a common isotropic value of permeability and porosity, equal to 1 × 10 -12 m 2 and 0.01, respectively.
The determination of the hydraulic properties for the CC complex is not straightforward. Indeed, the carbonate rock presents a series of fractures and discontinuities caused by tectonics and weathering processes [45], whose orientation and spatial distribution affects the hydraulic properties of the rock mass. Fractures induce anisotropy in the permeability tensor [46]. Local geometric and rheological anisotropies, joint density and orientation, and previous stress path significantly influence the hydraulic properties of rocks at the scale of tens of meters. Permeability and porosity from the literature show considerable variability, ranging from 10 -8 to 10 -14 m 2 and from 0.1 to 50%, respectively [47,48]. Model parameters depend on the size of the modelled domain. Therefore, since the modelling of each discrete fracture is unfeasible at the chosen simulation scale, an equivalent continuum approach is adopted [49]. It is assumed that each stratum is continuous, and we derive the equivalent hydraulic parameters from the properties of both the intact rock and joints. In particular, we followed a geomechanical approach [50,51] to characterise the equivalent rock mass permeability of the CC complex according to the following steps: (i) Identification of the fracturing state of the rock mass using geomechanical surveys (ii) Definition of a Representative Element Volume (REV) [37], typical of the rock-mass fracturing state and density (iii) Simulation of the fluid flow inside the discrete fractures of the REV and calculation of the equivalent permeability [51] (iv) Assignment of the estimated permeability to the CC complex of the large-scale numerical model in Figure 4 Geomechanical surveys allow analysing the fracturing state of a rock mass by identifying one or more families of rock joints, constituted by planes parallel to each other, i.e., the joint set. For hydromechanical purposes, the joint has been assumed smooth and planar. Thus, only the fracture's aperture and the spacing have been measured during the geomechanical surveys; the aperture is the distance between the walls of the single joint, while the spacing is the orthogonal distance between two joints of the same set. The measures involved the carbonate rock masses of the Eastern Aurunci Mts at four locations (coloured triangles in Figure 1(b)) and provided minimum and maximum values of aperture and spacing for every recognised joint set.

Geofluids
The performed geomechanical surveys ( Figure 5) highlight a well-organised fracturing setting with several joint sets that take origin from the stratigraphic and tectonic features of the investigated rock mass. Four joint sets have been defined for the investigated area, as shown in Table 2, with a diffuse high-angle dip between 70 and 90 degrees. The spacing and the aperture of these joint sets are variable and depend on the local heterogeneities that affect the rock mass, varying from 0.05 to 1.5 cm and from 5 to 120 cm, respectively. The rock also has bedding, with a subhorizontal dip and an average spacing and aperture of 30 cm and 0.05 cm, respectively.
According to the estimated fracturing state of the carbonate rock mass, we constructed a 2D REV ( Figure 6). The REV allows relating the microscopic hydraulic properties of the fractures to the macroscopic ones used in the equivalent continuum model of Figure 4.
The REV's x-and y-axes are oriented like those of the 2D numerical cross-section in Figure 4. Therefore, the x-axis corresponds to the horizontal fluid flow in the NW-SE direction (Figure 4), and the y-axis corresponds to the vertical fluid flow. The REV dimensions, i.e., 4 5 m × 4 5 m, are selected in order to represent the rock-mass fracture field adequately. Indeed, assuming a larger REV does not change the average fracture density inside the volume, as well as the REV's equivalent permeability. The REV's internal boundaries represent the discrete fractures, whose aperture and spacing are selected according to the results of the geomechanical surveys ( Table 2). The NW-SE joint set is not included since it is parallel to the modelled 2D cross-section ( Figure 4). The orientation and the spacing of the fractures have been assumed equal to the observed average values, while the aperture has been varied between the observed maximum and minimum values, respectively ( Table 2).   8 Geofluids The macroscopic equivalent permeability of the REV is then calculated by assuming the simultaneous presence of the bedding and the joint sets of Table 2 and simulating the flow numerically through the discrete fractures [37], according to the procedure proposed by Lancia et al. [51]. The fluid flow inside the fractures obeys to the Darcy law (equations (5) and (6)), while the flow through the porous matrix is neglected by assuming an arbitrarily low permeability. The permeability of the single fracture (κ f ) is calculated according to equation (9) and depends on the fracture's minimum and maximum aperture in Table 2. The flow through the rock fractures has been simulated alternatively along the x-and y-direction, respectively. Boundary conditions consist of an imposed pressure difference (p 1 and p 2 < p 1 ) along two opposite faces of the block, i.e., boundaries A and B for fluid flow along the x-direction and boundaries C and D for fluid flow along the y-direction ( Figure 6). The remaining boundaries are assumed impermeable. A porous-equivalent macroscopic permeability is then calculated at the end of the seepage flow analysis along the x (κ 0xx ) and y (κ 0yy ) directions by measuring the mean discharge along the outflowing boundary of the REV in steady-state conditions and applying the Darcy law (equation (3)), solved with respect to the permeability value.
The estimated equivalent permeability tensor is then applied to the CC complex of the 2D numerical model in Figure 4 to estimate the steady-state hydrothermal flow field. However, the macroscopic permeability tensor from the REV is representative of the near-surface fracture field only. At depth, the permeability of the Earth's crust generally decreases nonlinearly [48,52,53] because of the increase in lithostatic load and the resulting fracture closure. Thus, according to the literature [54][55][56], we assumed an exponential decrease of the permeability tensor of the CC complex with depth according to the following equation: In equation (11), z (m) is the height with respect to the mean sea level (Figure 4), κ 0ii (for i = 1, 2) is the permeability tensor calculated from the REV analysis, κ Lii <κ 0ii is the asymptotic permeability value at depth, α (-) is an exponential decay index, and β (m) is a threshold height, assumed equal to the average sea level, i.e., β = 0. The latter value accounts for a negligible permeability variation for the first 1000 meter depth inside the Eastern Aurunci Mts. Indeed the high fracture density and the presence of karst conduits affect the permeability more than the increasing stress [8]. Assuming higher β values does not change the results substantially.
The unknown parameters, i.e., the temperature T max (K) and gas concentration c 0 (mol/m 3 ) at the bottom of the model, the exponential decay coefficient α (-), the CC complex porosity ɛ p (-), the discrete fractures aperture (Figure 4)  9 Geofluids optimisation procedure [57]. The unknowns are varied between predefined ranges and selecting those values that minimise the sum of squared residuals (SSR) between the measured temperatures at selected water wells and springs and the corresponding modelled temperatures at the same positions: In equation (12), T o and T m are observed (i.e., measured) and modelled temperatures, respectively.

Equivalent Permeability of the Carbonate Complex.
The results of the numerical simulation of fluid flow inside the REV show an overall homogeneous steady-state flow through the discrete fractures along both the x and y directions (black arrows in Figures 7(a) and 7(b)). The macroscopic permeability of the REV is then estimated according to the Darcy law (equation (3)) as the ratio between the mean flow discharge, orthogonal to the flux direction, and the imposed pressure gradient. Since we assumed the maximum and minimum fracture aperture of Table 2, the final permeability tensor κ 0ii (13) is given by the mean value between the permeability tensors calculated assuming the minimum and maximum fracture aperture along the x-and y-direction, respectively. The estimated equivalent permeability tensor from the REV (13), representative of the near-surface fracture field, is then assigned to the CC complex of the 2D numerical model (Figure 4) through equation (11) to simulate the steady-state hydrothermal flow field.  Figure S1 in the supplementary material. The results of the numerical simulation are shown in Figure 8. The computed steady-state temperature distribution (Figure 8(a)) is maximum close to the heat source at the bottom right boundary of the model and then gradually decreases moving away from the imposed heat source. The estimated mean temperature inside the deep hydrostructure (the green dashed polygon in Figure 4) is approximately 169°C and agrees with the temperature interval of approximately 150°C -170°C estimated by geothermometers [8,10].
The comparison between the measured and the optimally computed temperatures is shown in Table 3. The water wells and springs selected for the optimisation procedure are those closest to the modelled A-B cross-section in Figure 1(b). There is a general agreement between the available temperature data and the modelled values. The measured temperature variability at depth for water wells W15 to W26 cannot be reproduced within our 2D model;   thus, we assumed a mean adiabatic value for the comparison with the modelled values. For the Gallo well, the computed temperature at the bottomhole reasonably agrees with the measured value at the same depth. The pore fluid diffusion (Figure 8(b)) develops entirely inside the CC complex, while the other lithologies, being less permeable, behave as a barrier for fluid flow. A dual fluid circulation develops inside the CC complex with average velocities of 1 × 10 -13 m/s, peaking up to 1 × 10 -9 m/s. The first one generates in the shallow aquifer of the Eastern Aurunci Mts (left side of the model in Figure 8(b)), and it is continuously fed by the imposed pressure gradient inside the karst aquifer. The other one develops in the deep aquifer below the Roccamonfina volcano (right side of the model in Figure 8(b)), and it is sustained by the convective loops triggered by the strong temperature gradient at depth. Because of this dual fluid circulation, the thermal and chemical features of the water flowing at the final delivery point, i.e., the Suio wells and springs (Figure 8(b)), are the result of the heat exchange between the cold water from the shallow karst circuit and the hot water from the deep aquifer.
Finally, the computed spatial distribution of the steadystate gas concentration and flux is shown in Figure 8(c). The deep reservoir water shows a maximum gas concentration up to 45 mol/m 3 because of the continuous feeding in diluted gases done by the nearby volcanic rocks, while the shallow reservoir on the left does not contain dissolved gases. The gas concentration gradually reduces to zero moving toward the ground surface. However, at the delivery point of the Suio wells and springs, the mixing of shallow water with the deep, gas-rich water produces a gas concentration of approximately 3-10 mol/m 3 .

Discussion
The simulation results unveiled the complexity of the Suio hydrothermal system. The applied diffusion-thermal approach identified the presence of a dual fluid flow circulation inside the CC complex. To the NW of the modelled section, the pore pressure gradient inside the Eastern Aurunci reservoir feeds the fluid flow of meteoric karst water. To the SE, the sharp temperature gradient sustains stationary convective loops, where the water increases its temperature and enriches in gas and non-karst-origin ions (Figures 8(b) and 8(c)). As a result, the water flowing at springs and intercepted by water wells along the Garigliano river (Figure 1(b)) shows temperatures, hydrochemical composition, and gas concentration that depend on the prevailing flow circulation system. At the NE and SW limits of the Suio hydrothermal basin (springs S2 and S12 to S19 in Figure 1(b)), the cold groundwater, with dominantly karst ions [8], originates from the surficial karst flow system, since the volcanic heat source is far and the convective loops are negligible. Conversely, moving towards the Roccamonfina volcano (springs S3 to S11 in Figure 1(b)), the measured increase in water temperature and non-karst-origin ions, together with the increase in dissolved gases testifies the prevailing action of deep convective loops.
The diffusion process is governed mainly by pressure and temperature gradients and to a lesser extent by the presence of dissolved gases. The latter is required to explain the apparent artesian behaviour of water wells located at the SE limit of the Eastern Aurunci hydrostructure (W1 to W27 in Figures 1(b) and 3(b)-3(d)), which is caused by the progressive degassing of diluted gas.
The modelled discontinuities dissecting the CC complex, i.e., the red segments in Figure 4, do not substantially modify the fluid flow at the scale of the performed analysis. The optimised fracture thickness associated with those discontinuities (d f = 2 63 × 10 −7 m) is unrealistically small, and the corresponding fracture permeability tensor does not modify the fluid diffusion substantially. Indeed, further models neglecting the discrete fractures(a and b in Figure 4) confirmed the results shown in Figure 8. Instead, the anisotropic equivalent permeability of the CC complex ( Figure S1) plays a fundamental role in controlling the fluid circulation, the temperature distribution, and the transport of diluted gases at the Suio wells and springs. Indeed, preliminary parametric analyses (not presented), assuming a constant isotropic permeability, provided a steady-state temperature distribution, fluid flow paths, and gas concentrations inconsistent with the hydrogeological observations [58].
The diffusion process occurs entirely in the CC complex, whose permeability governs fluid pressure, velocity, and thermal convection. Assigning a different permeability to the VC, FC, GC, and MC complexes (Figure 4) does not modify the fluid diffusion significantly. Because permeability actively  12 Geofluids controls the fluid diffusion in porous media, an accurate estimation of the fracture-related permeability tensor in karst structures is required for a better prediction of the behaviour of complex hydrothermal systems. Despite the overall agreement between the modelling results and the conceptual hydrogeological model of the area [8], it is worth discussing the significance of the results regarding the modelling assumptions. In our modelling, the permeability tensor of the CC complex has been estimated using geomechanical surveys. Such an approach reconstructs the rock fracture heterogeneities by local measurements that often are unevenly distributed over the study area, because vegetation or detritus may hinder the carbonate unit. Moreover, without explorable karst conduits and caves, all the measures are performed along the topographic surface. Thus, the obtained joint sets could not be representative of the fracturing condition in depth. For this reason, we assumed an exponential decay of permeability with depth (equation (11)). Though the assumed decay is difficult to validate with in situ data, it is consistent with the available literature [48,52].
The fracture permeability in both the 2D numerical model and the REV has been estimated according to the cubic law [34] (equation (9)), which is valid under laminar fluid flow only. Laminar conditions are not acceptable in case of flux through extensive fractures or karst conduits where the flow is locally turbulent. Given the regional character of the work and considering that karst processes are young and not well-developed in the studied area [8], turbulent conditions have not been assumed in our analysis.
The negligible effect of the modelled discrete fractures on the overall fluid circulation is related to the scale of the performed modelling. At the macroscale, the equivalent fracture permeability of the rock mass affects the observed temperature distribution and the dual fluid flow, rather than the presence of discrete discontinuities. However, at medium-to-small scale, the presence of conduits could influence the fluid flow and the heat transfer. In case of local scale analyses, specific laboratory and field tests are required to investigate the relationship between fractures and fluid circulation and the transition from laminar to turbulent flow.
The assumed 2D geometry does not allow estimating water temperatures, discharge, and gas concentrations at each of the springs and water wells along the Garigliano river (Figure 1(b)). However, the performed modelling focuses more about the process characterisation than effective resource evaluation [59]. Therefore, the general patterns of temperatures, fluid flow, and gas concentration (Figure 8) have been considered for describing and validating the investigated scenario, rather than absolute magnitudes. A 3D model is necessary for a proper spatial estimation of each quantity and correct well and spring exploitation.
Overall, the results of the performed numerical model confirm the hydrogeological setting proposed by Saroli et al. [8] and agree with the available literature over the study area. Indeed, the estimated gas concentration at the delivery point of Suio wells and springs (3-10 mol/m 3 ) is in agreement with the approximately 4-6 mol/m 3 measured at some springs along the Garigliano river [60]. Moreover, the computed temperature distribution at depth approximately follows the estimated subsurface temperature distributions for the whole Italian territory at 1, 2, and 3 kilometers below the sea level [2] ( Figure S2 in the supplementary material).
The findings of the numerical model confirm the presence of a medium-to-low temperature reservoir. By this model, preliminary quantitative evaluation of the groundwater resource and its geothermal potential can be derived with minimum future efforts. Indeed, the investigated hydrothermal system represents a significant energy resource, potentially exploitable for both urban heating systems and industrial processes. In particular, this model is suitable to predict temperature variations of springs and wells if an increasing water withdrawal is performed, or in a different climate scenario where the recharge from carbonate aquifer is shortened. More generally, the calibrated numerical model can be used as a predictive tool for the simulation of potential exploitation scenarios and sustainable use of the geothermal resource.

Conclusions
In this work, the complex behaviour of the Suio hydrothermal system has been validated employing a hybrid continuum-discrete numerical model, calibrated with available and new data.
Simulations coupled fluid diffusion, thermal convection, and transport of diluted species in porous media, as well as the flow through discrete fractures. The simulation results show that the Suio hydrothermal activity, i.e., the temperature, chemical composition, and diluted gas of the hot springs and wells, is linked to (i) a heat source at depth, related to the peri-Tyrrhenian volcanism, responsible for the strong temperature gradient (ii) the mutual interaction between the dual fluid flow circulation inside the CC complex; one sustained by the hydraulic gradient inside the Eastern Aurunci Mts and the other feed by convective loops (iii) the fracture-related anisotropic permeability of the CC complex, which regulates both the heat transfer, fluid diffusion, and transport of diluted gas Such results corroborate and numerically quantify the hydrogeological conceptual model of the Suio hydrothermal area [8], allowing for an easy future first-order evaluation of sustainable exploitation and planning of the geothermal resource.

Data Availability
The data used to support the findings of this study are included within the article.