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

Saroli, Michele; Albano, Matteo; Giovinco, Gaspare; Casale, Anna; Dell’Isola, Marco; Lancia, Michele; Petitta, Marco

doi:https://doi.org/10.1155/2019/5485068

Geofluids

On this page

Abstract Introduction Data and Methods Results Discussion Conclusions Data Availability Conflicts of Interest Acknowledgments Supplementary Materials References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 5485068 | https://doi.org/10.1155/2019/5485068

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

Michele Saroli,^1,2Matteo Albano ,²Gaspare Giovinco,¹Anna Casale,³Marco Dell’Isola,¹Michele Lancia,^4,5and Marco Petitta⁶

Guest Editor: Philippe Pasquier

Received23 Jan 2019

Accepted21 Mar 2019

Published16 May 2019

Abstract

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.

1. Introduction

In recent decades, the worldwide growth of energy demand and the increase in CO₂ 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.

(a)

(b)

(c)

Figure 1

(a) Hydrogeological map of the Eastern Aurunci Mts and Roccamonfina volcano modified from [8]. The inset shows the mean temperature at a 3000-meter depth [2]. (b) Detail of the Suio Terme area (the back dashed rectangle in (a)). The numbers refer to the surveyed springs (S) and water wells (W). (c) Geological and hydrogeological section of the Suio area, passing through the A-B cross-section line in (a). Key to the legend: (1) alluvial complex (AC); (2) volcanic complex (VC); (3) sandy-conglomeratic complex (GC); (4) turbidite complex (FC); (5) carbonate complex (CC), the light-blue circles indicate the presence of gas; (6) lower crust or metamorphic rocks (MC); (7) normal faults: (A) certain, (B) inferred; (8) thrust faults: (A) certain, (B) inferred; (9) inferred faults ((c) only): (A) normal and (B) thrust; (10) Caldera nest ((a) only); (11) main springs ((a) only); (12) springs with temperature (A) °C, (B) °C, (C) °C ((b and c) only); (13) linear springs; (14) water wells with temperature (A) °C, (B) °C, (C) °C, (D) °C; (15) groundwater flux (blue arrow) and piezometric line (dashed blue curves). The numbers indicate the groundwater elevation in meters a.s.l. (a and b); (16) groundwater table ((c) only); (17) cold (A) and thermal (B) groundwater fluxes ((c) only); (18) the position of the Gallo well (GW) 85-1 [7] ((a and c) only); and (19) location of the geomechanical stations ((b) only).

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–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–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.

2. 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 high-discharge 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 the Suio hydrothermal basin (S2 and S19 in Figure 1(b)), the groundwater is cold (°C) with a predominance of dissolved karst ions (Ca²⁺, HCO₃^-). 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₄^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₂ and subordinately H₂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].

2.2. 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³) is the density of the fluid, (J/kg K) is the fluid heat capacity at constant pressure, (J/m³ K) is the effective volumetric heat capacity at constant pressure, (K) is the temperature, (m/s) is the fluid velocity field, (W/m³) is the heat source, and (W/m K) is the effective thermal conductivity of the solid-fluid system. is given by the weighted arithmetic mean of fluid and porous matrix conductivities: where and are the thermal conductivity of the solid and fluid, respectively, and (-) is the volume fraction of the solid, given by , where (-) 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 (m/s) depends on the permeability (m²), the fluid’s dynamic viscosity (Pa s), and is driven by a pressure gradient (Pa). Darcy’s law is then combined with the continuity equation: where (kg/m³) is the fluid density and (m³) 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 denotes the gradient operator restricted to the fracture’s tangential plane, (m) is the fracture thickness, and (m²) and (-) 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: where (mol/m³) denotes the concentration of species in the fluid, (mol/m³ s) is a reaction rate expression for the species, (m/s) is the solvent velocity field, and (m²/s) is the molecular diffusion coefficient of the diluted species, assumed temperature-dependent according to the following equation: where m²/s is the average molecular diffusion in water, and (K) is the temperature [26–28].

2.3. Finite Element Model Setup

The conceptual cross section of Figure 1(c) [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].

Figure 4

A 2D finite element model derived from the hydrogeological section in Figure 1(c). The model shows the location of the temperature measurement points of the Suio springs and at the bottom of the Gallo well and the Suio water wells along the Garigliano river. The area bordered by the green dashed polygon represents the deep reservoir. The dashed blue line represents the piezometric level inside the carbonates. The two red segments (a and b) represent the modelled discrete fractures. For the nomenclature of each layer, please refer to the caption of Figure 1.

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 [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 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 W/m² 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 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 mol/m³.

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 (equations (5) and (6)) and the relation between fracture aperture and its corresponding permeability (m²) in the direction parallel to the fluid flow is derived from the well-known cubic law [34]: where (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 (m²) is the permeability tensor with respect to the global coordinate system and is the positive counter-clockwise rotation angle between the global -axis of Figure 4 and the modelled discrete fractures.

2.4. 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₂ 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–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 ^-12 m² 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² 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. 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 - and -axes are oriented like those of the 2D numerical cross-section in Figure 4. Therefore, the -axis corresponds to the horizontal fluid flow in the NW-SE direction (Figure 4), and the -axis corresponds to the vertical fluid flow. The REV dimensions, i.e., , 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).

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 () 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 - and -direction, respectively. Boundary conditions consist of an imposed pressure difference ( and ) along two opposite faces of the block, i.e., boundaries A and B for fluid flow along the -direction and boundaries C and D for fluid flow along the -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 () and () 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–56], we assumed an exponential decrease of the permeability tensor of the CC complex with depth according to the following equation:

In equation (11), (m) is the height with respect to the mean sea level (Figure 4), (for ) is the permeability tensor calculated from the REV analysis, 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., . 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 (K) and gas concentration (mol/m³) at the bottom of the model, the exponential decay coefficient (-), the CC complex porosity (-), the discrete fractures aperture (Figure 4) (m), and the asymptotic permeability values (m²) (for ) are estimated with the Nelder-Mead 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), and are observed (i.e., measured) and modelled temperatures, respectively.

3. Results

3.1. 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 and 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 (13) is given by the mean value between the permeability tensors calculated assuming the minimum and maximum fracture aperture along the - and -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.

(a)

(b)

3.2. Steady-State Model

The estimated parameters at the end of the optimization procedure are °C, mol/m³, , , m, and m², respectively. The horizontal and vertical permeability distributions inside the CC complex are shown in 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].

(a)

(b)

(c)

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 ^-13 m/s, peaking up to ^-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 steady-state gas concentration and flux is shown in Figure 8(c). The deep reservoir water shows a maximum gas concentration up to 45 mol/m³ 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³.

4. 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 ( 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 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³) is in agreement with the approximately 4-6 mol/m³ 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.

5. 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.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

The authors would like to thank Dr Pierino De Felice (INMRI-ENEA–Agenzia nazionale per le tecnologie, l’energia e lo sviluppo economico sostenibile) who partially financed this research work and Dr Geol. Francesco Nolasco who provided the pictures in Figures 3(c) and 3(d). Perceptually uniform colourmaps are used in specific figures [61].

Supplementary Materials

The supplementary material contains two figures showing the permeability distribution inside the 2D numerical model and the comparison between the modelled temperatures of Figure 8(a) and the values estimated at 1000, 2000, and 3000 meters below the sea level by Cataldi et al. [2]. (Supplementary Materials)

References

M. Procesi, B. Cantucci, M. Buttinelli, G. Armezzani, F. Quattrocchi, and E. Boschi, “Strategic use of the underground in an energy mix plan: synergies among CO₂, CH₄ geological storage and geothermal energy. Latium region case study (central Italy),” Applied Energy, vol. 110, pp. 104–131, 2013.
View at: Publisher Site | Google Scholar
R. Cataldi, F. Mongelli, P. Squarci, L. Taffi, G. Zito, and C. Calore, “Geothermal ranking of Italian territory,” Geothermics, vol. 24, no. 1, pp. 115–129, 1995.
View at: Publisher Site | Google Scholar
B. Billi, G. Cappetti, and F. Luccioli, “Enel activity in the research, exploration and exploitation of geothermal energy in Italy,” Geothermics, vol. 15, no. 5-6, pp. 765–779, 1986.
View at: Publisher Site | Google Scholar
P. Muffler and R. Cataldi, “Methods for regional assessment of geothermal resources,” Geothermics, vol. 7, no. 2-4, pp. 53–89, 1978.
View at: Publisher Site | Google Scholar
D. Cinti, M. Procesi, and P. Poncia, “Evaluation of the theoretical geothermal potential of inferred geothermal reservoirs within the Vicano–cimino and the Sabatini volcanic districts (central Italy) by the application of the volume method,” Energies, vol. 11, no. 1, 2018.
View at: Publisher Site | Google Scholar
G. Buonasorte, R. Cataldi, T. Franci et al., Previsioni di crescita della geotermia in Italia fino al 2030–Per un Nuovo Manifesto della Geotermia Italiana, Pacini, UGI Special Publication, 2011.
M. D. Watts, “Geothermal exploration of Roccamonfina volcano, Italy,” Geothermics, vol. 16, no. 5-6, pp. 517–528, 1987.
View at: Publisher Site | Google Scholar
M. Saroli, M. Lancia, M. Albano et al., “A hydrogeological conceptual model of the Suio hydrothermal area (central Italy),” Hydrogeology Journal, vol. 25, no. 6, pp. 1811–1832, 2017.
View at: Publisher Site | Google Scholar
C. Nunziata and F. Gerecitano, “VS crustal models of the Roccamonfina volcano and relationships with Neapolitan volcanoes (southern Italy),” International Journal of Earth Sciences, vol. 101, no. 5, pp. 1371–1383, 2012.
View at: Publisher Site | Google Scholar
F. D’Amore, A. Didomenicoantonio, and S. Lombardi, “Considerazioni geochimiche e geotermometriche sul sistema idrotermale di Suio (Campania),” Geologica Romana, vol. 31, pp. 319–328, 1995.
View at: Google Scholar
A. Taillefer, R. Soliva, L. Guillou-Frottier, E. Le Goff, G. Martin, and M. Seranne, “Fault-related controls on upward hydrothermal flow: an integrated geological study of the Têt fault system, Eastern Pyrénées (France),” Geofluids, vol. 2017, Article ID 8190109, 19 pages, 2017.
View at: Publisher Site | Google Scholar
R. Therrien and E. A. Sudicky, “Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media,” Journal of Contaminant Hydrology, vol. 23, no. 1-2, pp. 1–44, 1996.
View at: Publisher Site | Google Scholar
A. M. Shapiro and J. Andersson, “Steady state fluid response in fractured rock: a boundary element solution for a coupled, discrete fracture continuum model,” Water Resources Research, vol. 19, no. 4, pp. 959–969, 1983.
View at: Publisher Site | Google Scholar
T. C. Rasmussen and D. D. Evans, Fluid Flow and Solute Transport Modeling through Three-Dimensional Networks of Variably Saturated Discrete Fractures, Nuclear Regulatory Commission, 1989.
J. C. S. Long, J. S. Remer, C. R. Wilson, and P. A. Witherspoon, “Porous media equivalents for networks of discontinuous fractures,” Water Resources Research, vol. 18, no. 3, pp. 645–658, 1982.
View at: Publisher Site | Google Scholar
B. Faybishenko, P. A. Witherspoon, and S. M. Benson, Dynamics of Fluids in Fractured Rock; Geophysical Monograph Series, American Geophysical Union, Washington, DC, USA, 2000.
S. P. Neuman, “Trends, prospects and challenges in quantifying flow and transport through fractured rocks,” Hydrogeology Journal, vol. 13, no. 1, pp. 124–147, 2005.
View at: Publisher Site | Google Scholar
P. Lachassagne, S. Ahmed, C. Golaz et al., “A methodology for the mathematical modelling of hard-rock aquifers at catchment scale based on the geological structure and the hydrogeological functioning of the aquifer,” in XXXI IAH Congress, Munich, K.-P. Seiler and S. Wohnlich, Eds., pp. 367–370, New Approaches Characterising Groundwater Flow, AA Balkema Publishers, 2001.
View at: Google Scholar
L. Jing and J. A. Hudson, “Numerical methods in rock mechanics,” International Journal of Rock Mechanics and Mining Sciences, vol. 39, no. 4, pp. 409–427, 2002.
View at: Publisher Site | Google Scholar
D. Cosentino, I. Federici, P. Cipollari, and E. Gliozzi, “Environments and tectonic instability in central Italy (Garigliano Basin) during the late Messinian Lago–Mare episode: new data from the onshore Mondragone 1 well,” Sedimentary Geology, vol. 188-189, pp. 297–317, 2006.
View at: Publisher Site | Google Scholar
F. Ippolito, F. Ortolani, and M. Russo, “Struttura marginale tirrenica dell’Appennino campano: reinterpretazione di dati di antiche ricerche di idrocarburi,” Memorie della Società Geologica Italiana, vol. 12, pp. 227–250, 1973.
View at: Google Scholar
M. Bernabini, D. Bucci, L. Orlando, M. Parotto, M. M. Tiberti, and M. Tozzi, “Nuovi vincoli gravimetrici 3D in Italia centro-meridionale per un modello della crosta profonda,” Memorie Società Geologica Italiana, vol. 57, no. 1, pp. 335–342, 2002.
View at: Google Scholar
C. Savelli, “Late Oligocene to recent episodes of magmatism in and around the Tyrrhenian Sea: implications for the processes of opening in a young inter-arc basin of intra-orogenic (Mediterranean) type,” Tectonophysics, vol. 146, no. 1-4, pp. 163–181, 1988.
View at: Publisher Site | Google Scholar
E. Cuoco, G. Verrengia, S. de Francesco, and D. Tedesco, “Hydrogeochemistry of Roccamonfina volcano (southern Italy),” Environment and Earth Science, vol. 61, no. 3, pp. 525–538, 2010.
View at: Publisher Site | Google Scholar
COMSOL Inc, COMSOL Multiphysics® Reference Manual V.5.4, COMSOL AB, Stockholm, Sweden, 2018.
J. Welty, G. L. Rorrer, and D. G. Foster, Fundamentals of Momentum, Heat and Mass Transfer, John Wiley & Sons Inc, 2014.
F. Zhao, T. J. Armstrong, and A. V. Virkar, “Measurement of O₂-N₂ effective diffusivity in porous media at high temperatures using an electrochemical cell,” Journal of the Electrochemical Society, vol. 150, no. 3, p. A249, 2003.
View at: Publisher Site | Google Scholar
S. P. Cadogan, G. C. Maitland, and J. P. M. Trusler, “Diffusion coefficients of CO₂ and N₂ in water at temperatures between 298.15 K and 423.15 K at pressures up to 45 MPa,” Journal of Chemical & Engineering Data, vol. 59, no. 2, pp. 519–525, 2014.
View at: Publisher Site | Google Scholar
F. Mostardini and S. Merlini, “Appennino centro meridionale: sezioni geologiche e proposta di modello strutturale,” Memorie della Società Geologica Italiana, vol. 35, pp. 177–202, 1986.
View at: Google Scholar
P. Celico, “Idrogeologia dei massicci carbonatici, delle piane quaternarie e delle aree vulcaniche dell’Italia centromeridionale (Marche, Lazio meridionale, Abruzzo, Molise e Campania),” Cassa per il Mezzogiorno, vol. 4, pp. 1–203, 1983.
View at: Google Scholar
C. F. Boni, P. Bono, and G. Capelli, “Schema idrogeologico dell’Italia centrale,” Memorie della Società Geologica Italiana, vol. 35, pp. 991–1012, 1986.
View at: Google Scholar
T. F. M. Chui and D. L. Freyberg, “Implementing hydrologic boundary conditions in a multiphysics model,” Journal of Hydrologic Engineering, vol. 14, no. 12, pp. 1374–1377, 2009.
View at: Publisher Site | Google Scholar
T. L. Bergman, A. S. Lavine, F. P. Incropera, and D. P. DeWitt, Fundamentals of Heat and Mass Transfer, John Wiley & Sons Inc, 2011.
P. A. Witherspoon, J. S. Y. Wang, K. Iwai, and J. E. Gale, “Validity of cubic law for fluid flow in a deformable rock fracture,” Water Resources Research, vol. 16, no. 6, pp. 1016–1024, 1980.
View at: Publisher Site | Google Scholar
J. Huang and A. Ghassemi, “A poroelastic model for evolution of fractured reservoirs during gas production,” Journal of Petroleum Science and Engineering, vol. 135, pp. 626–644, 2015.
View at: Publisher Site | Google Scholar
B. Hölting and W. G. Coldewey, “Hydrogeology,” in Springer Textbooks in Earth Sciences, Geography and Environment, pp. 1–357, Springer Berlin Heidelberg, Berlin, Heidelberg, 2019.
View at: Publisher Site | Google Scholar
J. Bear, Dynamics of Fluids in Porous Media, American Elsevier Pub. Co., New York, NY, USA, 1972.
P. Froldi, Impianti Geotermici. Progettazione, Realizzazione e Controllo, Maggioli Editore, 2013.
G. Buonanno, A. Carotenuto, G. Giovinco, and N. Massarotti, “Experimental and theoretical modeling of the effective thermal conductivity of rough steel spheroid packed beds,” Journal of Heat Transfer, vol. 125, no. 4, p. 693, 2003.
View at: Publisher Site | Google Scholar
M. Dell’Isola, F. R. d’Ambrosio Alfano, G. Giovinco, and E. Ianniello, “Experimental analysis of thermal conductivity for building materials depending on moisture content,” International Journal of Thermophysics, vol. 33, no. 8-9, pp. 1674–1685, 2012.
View at: Publisher Site | Google Scholar
A. Frattolillo, G. Giovinco, M. C. Mascolo, and A. Vitale, “Effects of hydrophobic treatment on thermophysical properties of lightweight mortars,” Experimental Thermal and Fluid Science, vol. 30, no. 1, pp. 27–35, 2005.
View at: Publisher Site | Google Scholar
P. V. Sharma, Environmental and Engineering Geophysics, Cambridge University Press, Cambridge, UK, 1997.
View at: Publisher Site
S. Viaroli, L. Mastrorillo, F. Lotti, V. Paolucci, and R. Mazza, “The groundwater budget: a tool for preliminary estimation of the hydraulic connection between neighboring aquifers,” Journal of Hydrology, vol. 556, pp. 72–86, 2018.
View at: Publisher Site | Google Scholar
S. Viaroli, L. Mastrorillo, and R. Mazza, “Contribution of the Roccamonfina caldera to the basal volcanic aquifer recharge: first considerations,” Rendiconti Online della Società Geologica Italiana, vol. 41, pp. 95–98, 2016.
View at: Publisher Site | Google Scholar
P. Lachassagne, R. Wyns, and B. Dewandel, “The fracture permeability of hard rock aquifers is due neither to tectonics, nor to unloading, but to weathering processes,” Terra Nova, vol. 23, no. 3, pp. 145–161, 2011.
View at: Publisher Site | Google Scholar
D. T. Snow, “Anisotropie permeability of fractured media,” Water Resources Research, vol. 5, no. 6, pp. 1273–1289, 1969.
View at: Publisher Site | Google Scholar
F. Agosta, M. Prasad, and A. Aydin, “Physical properties of carbonate fault rocks, Fucino basin (central Italy): implications for fault seal in platform carbonates,” Geofluids, vol. 7, no. 1, 32 pages, 2007.
View at: Publisher Site | Google Scholar
T. Gleeson and S. Ingebritsen, “Crustal Permeability,” Tech. Rep., John Wiley & Sons, Ltd., Chichester, UK, 2012.
View at: Google Scholar
T. G. Sitharam, J. Sridevi, and N. Shimizu, “Practical equivalent continuum characterization of jointed rock masses,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 3, pp. 437–448, 2001.
View at: Publisher Site | Google Scholar
W. E. Bamford, C. M. Barton, B. MacMahon et al., “International society for rock mechanics commission on standardization of laboratory and field tests: suggested methods for the quantitative description of discontinuities in rock masses,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 15, no. 6, pp. 319–368, 1978.
View at: Publisher Site | Google Scholar
M. Lancia, M. Saroli, and M. Petitta, “A double scale methodology to investigate flow in karst fractured media via numerical analysis: the Cassino plain case study (central Apennine, Italy),” Geofluids, vol. 2018, Article ID 2937105, 12 pages, 2018.
View at: Publisher Site | Google Scholar
T. Gleeson, N. Moosdorf, J. Hartmann, and L. P. H. van Beek, “A glimpse beneath earth’s surface: GLobal HYdrogeology MaPS (GLHYMPS) of permeability and porosity,” Geophysical Research Letters, vol. 41, no. 11, pp. 3891–3898, 2014.
View at: Publisher Site | Google Scholar
X. Kuang and J. J. Jiao, “An integrated permeability-depth model for Earth’s crust,” Geophysical Research Letters, vol. 41, no. 21, pp. 7539–7545, 2014.
View at: Publisher Site | Google Scholar
X.-W. Jiang, X.-S. Wang, and L. Wan, “Semi-empirical equations for the systematic decrease in permeability with depth in porous and fractured media,” Hydrogeology Journal, vol. 18, no. 4, pp. 839–850, 2010.
View at: Publisher Site | Google Scholar
M. Albano, S. Barba, M. Saroli et al., “Aftershock rate and pore fluid diffusion: insights from the Amatrice-Visso-Norcia (Italy) 2016 seismic sequence,” Journal of Geophysical Research: Solid Earth, vol. 124, no. 1, pp. 995–1015, 2019.
View at: Publisher Site | Google Scholar
J. Hommel, E. Coltman, and H. Class, “Porosity–permeability relations for evolving pore space: a review with a focus on (bio-)geochemically altered porous media,” Transport in Porous Media, vol. 124, no. 2, pp. 589–629, 2018.
View at: Publisher Site | Google Scholar
J. A. Nelder and R. Mead, “A simplex method for function minimization,” The Computer Journal, vol. 7, no. 4, pp. 308–313, 1965.
View at: Publisher Site | Google Scholar
A. Casale, M. Saroli, G. Giovinco et al., “Hydrogeological study and numerical model of the Suio-Castelforte hydrothermal area (central Italy),” Rendiconti Online della Società Geologica Italiana, vol. 41, pp. 276–279, 2016.
View at: Publisher Site | Google Scholar
S. L. Barbour and J. Krahn, “Numerical modelling–prediction or process,” Geotechnical News, vol. 22, no. 4, pp. 44–52, 2004.
View at: Google Scholar
F. Frondini, C. Cardellini, S. Caliro, G. Beddini, A. Rosiello, and G. Chiodini, “Measuring and interpreting CO₂ fluxes at regional scale: the case of the Apennines, Italy,” Journal of the Geological Society, vol. 176, no. 2, pp. 408–416, 2019.
View at: Publisher Site | Google Scholar
F. Crameri, Scientific colour-maps, Zenodo, 2018.
View at: Publisher Site

Copyright

Copyright © 2019 Michele Saroli et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

957

Downloads

1128

Citations