Low Enthalpy Geothermal Systems in Structural Controlled Areas: A Sustainability Analysis of Geothermal Resource for Heating Plant (The Mondragone Case in Southern Appennines, Italy)

Abstract

In this study, the sustainability of low-temperature geothermal field exploitation in a carbonate reservoir near Mondragone (CE), Southern Italy, is analyzed. The Mondragone geothermal field has been extensively studied through the research project VIGOR (Valutazione del potenzIale Geotermico delle RegiOni della convergenza). From seismic, geo-electric, hydro-chemical and groundwater data, obtained through the experimental campaigns carried out, physiochemical features of the aquifers and characteristics of the reservoir have been determined. Within this project, a well-doublet open-loop district heating plant has been designed to feed two public schools in Mondragone town. The sustainability of this geothermal application is analyzed in this study. A new exploration well (about 300 m deep) is considered to obtain further stratigraphic and structural information about the reservoir. Using the derived hydrogeological model of the area, a numerical analysis of geothermal exploitation was carried out to assess the thermal perturbation of the reservoir and the sustainability of its exploitation. The effect of extraction and reinjection of fluids on the reservoir was evaluated for 60 years of the plant activity. The results are fundamental to develop a sustainable geothermal heat plant and represent a real case study for the exploitation of similar carbonate reservoir geothermal resources.

1. Introduction

The total (electric and thermal) potential energy of a geothermal reservoir depends on many variables such as thermal regime, geological structure, and hydrogeological properties. Recently, the use of groundwater, extracted via open-loop systems, is increasingly being considered as an efficient means for heating and cooling of buildings, especially in the heat pump systems.

This research concerns the geothermal carbonate reservoir, which makes the bedrock in the alluvial Mondragone plain, Southern Italy. The site has been extensively studied through the VIGOR project (coordinated by the Italian National Council Research and sponsored by the Italian Ministry of Economic Development) dedicated to the evaluation of geothermal potential in four Italian Regions: Puglia, Calabria, Campania, and Sicily. In this plain, a well-doublet open-loop plant has been designed for the district heating of seven public schools of Mondragone town [1].

Low-temperature geothermal fields have been exploited over decades for industrial applications in Iceland, Hungary, China, Turkey, France, Germany, Russia, and other countries [2], and in recent times many authors propose the use of geothermal energy standalone or coupled with other renewable energy sources [3,4,5]. As an example, a classical low-temperature geothermal field of meteoric origin, in Kamchatka, Russia, has extracted thermal water since 1966 mainly in the mode of artesian flow to supply numerous swimming pools, district heating of two villages, greenhouse farming, and fish breeding [6]. The long-term exploitation (20 years) has been assessed finding an insignificant temperature decrease in the geothermal reservoir (0.5 °C) [2]. The geothermal power potential in a hot spring close to the Municipality of Isa, in Japan, has been analyzed through the gravity data estimating a power density of 30.4 kW_e/km² for a 20 year exploitation [7]. Another interesting investigation has been carried out in a volcanic research site in the North German Basin, where a well-doublet with one injection and one production well has been explored for future geothermal energy production [8]. Alternating the injection of low and high flow rates of water a rapid increase in the water level of an adjacent well has been observed as well as an increase of the overall productivity of the treated well [9].

With the exception of active volcanic areas, carbonate aquifers are the most important geothermal reservoirs in which low-to-medium temperatures are usually reached at significant depths from the ground level (>1000 m) (e.g., [10,11,12,13,14,15,16]).

In many countries (i.e., Southern Italy), fault-controlled carbonate extensional domains may present temperatures in the range of low or middle enthalpy (T < 150 °C), at relatively shallow depths (<0.5 km) [14,17,18,19,20,21]. In fact, fault-controlled geothermal sites are accompanied by hydrothermal activity, hot fluid circulation, and mineralization processes [22,23]. Fracturing and chemical dissolution of carbonates increase the permeability of rocks, enhancing the advection of hot fluids and generating heat and mass transport [24,25,26,27,28]. The setting of carbonate bedrock in Mondragone plain corresponds to this hydrogeological model.

In recent times, the installation of innovative pilot geothermal electric and heating power plants has been incentivized in Italy. In this framework, prior to drilling activities and plant design, knowing scientific-technical data related to the potential of geothermal reservoirs and the sustainability of the utilization represents a crucial task to assess the economic feasibility of the geothermal exploitation and the development of future project plants [29,30].

Other information can be obtained in such a way, such as the relevance of system boundary conditions during long-term utilization, the interference of fluids extraction and reinjection during production, and the effectiveness of geothermal production systems. Numerical simulation is recognized as a fundamental tool for the elaboration and assessment of using geothermal energy [31], which is otherwise considered to be highly risky [32,33,34].

In recent years, different multiphase simulation tools have been widely used to model the effects induced by the exploitation of geothermal resources [35,36,37]. Among them, the code TOUGH2® (Transport Of Unsaturated Groundwater and Heat) has been firstly applied to Wairakei (New Zealand) geothermal field [38,39,40] and subsequently to several other geothermal fields (i.e., [41,42]). Additionally, the Los Alamos National Laboratory Finite Element Heat and Mass Transfer (FEHM®) code has a good record of numerical modeling studies, mainly related to enhanced geothermal systems (EGS) [43,44,45,46], together with the U.S. Geological Survey code HYDROTHERM® [47,48,49]. The commercial software COMSOL®, in particular, is a widely used simulator that has been adapted (as for example with the link to MATLAB®) [50] to applications including geothermal studies and, in several cases, the evaluation of fault influence [23,51,52]. As the heat and mass transfer in the reservoir depends on the internal properties of the reservoir and of the surrounding rocks, there has been an increased effort to model geothermal porous reservoirs. However up to now, few sustainability exploitation analyses on fractured carbonate rocks have been done [53,54] and they mainly focus on middle-high enthalpy and deep geothermal reservoirs, establishing several dynamic relations between the properties of the equivalent porous medium and fracture aperture. The exploitation of a low enthalpy geothermal system has been proposed by Cherubini et al. [55], who simulated a fractured geological system, composed of a single homogeneous layer with an inclined fault. They found that the width and permeability of the fault significantly affect fluid flow and thermal field, with a concentration of thermal perturbation within the fault plane.

2. Research Object Description

In the present paper, a sustainability analysis of the exploitation, by means of a well-doublet open-loop plant, of a low-temperature reservoir consisting of tectonically fractured and karstic limestone is carried out. This aquifer hosts an artesian groundwater body with a pressure larger than 1.70 bars, a wellhead temperature of 14 °C, and a natural water flow of 23.5 L/s. In the present work, two of the seven schools analyzed within the VIGOR Project are considered to be supplied. The characteristics of the system proposed to feed the two public schools are reported in

Table 1. The proposed system for the feeding of two public schools.

Thermal Power Plant	Network	Users
- Drawn water flow rates 6 L/s - One Titanium heat exchanger - Two Geothermal heat pumps with Tevap. = 25.0 °C - Total thermal power 160 kW - COP 4.00	- Principal Loop DN 100; branches DN 80 - Water flow rates 4.00 L/s - Network flow temperature 55.0 °C	- School Name “II Circolo” - School Name “Leonardo da Vinci” - Gym of School Name “Leonardo da Vinci”

.

A water flow of 6.00 out of 23.5 L/s was computed as sufficient to heat the considered schools (about 160 kW). The remainder is supposed to be by-passed and then mixed again to reinject the entire flow rate through a reinjection well [56]. A scheme describing the operation of the system is shown in Figure 1.

The sustainability of such geothermal exploitation has been analyzed by a numerical simulation, comprising a one fault case and two faults case, using the commercial software COMSOL Multiphysics® to assess the reservoir thermal perturbation (for 60 years of system operation for heating needs) due to extraction and reinjection of fluids at a rate suitable to produce the required thermal power.

Due to low enthalpy values and high water pressure effects on the reservoir, mechanical deformations have been considered negligible. A further element of interest in the modeling is due to the position of the wells (of the extraction/injection) in the plant. In fact, for logistical constraints, the injection well is not downstream of the extraction well, but almost flanked by this with respect to the flow direction of the groundwater body. The numerical results obtained allow development of a solid conceptual model for the Mondragone geothermal reservoir exploitation and open the possibility to future applications to similar geothermal reservoirs which are widespread in the Mediterranean region.

3. Hydrogeological Setting and Conceptual Model of the Geothermal Area

The area of interest is located on the coastal plain of Mondragone town at the bottom of Mt. Petrino carbonate hill (Figure 2).

This is at SW of Mt. Massico, which is a carbonate monocline trending SW bounded by high-angle normal faults along the margins of the Garigliano (NW), Volturno (SE), and Mondragone (SW) plains. These faults, several kilometers deep, have a throw estimated to be hundreds of meters [57,58,59]. Recent geophysical data have provided information about the bedrock of Mondragone plain. It is made by limestone and marly-arenaceous rocks and affected by several faults. The most important fault of this area crosses the whole plain from NW to SE and it is here coupled to other sub-parallel faults. The carbonate bedrock near Mt. Petrino has been found at a depth of –30 m below ground level (b.g.l.), under pyroclastic-alluvial deposits. From this zone, moving towards SE, the bedrock shows instead a regular depth towards the sea [1] (Figure 2).

Boreholes’ data in the plain have shown the presence of a powerful bank of Campanian Ignimbrite tuff almost continuous across the plain. This tuff has low permeability and therefore separates the pyroclastic-alluvial aquifer (very poor) above the tuff, from the underlain alluvial deposits, deeper and more permeable, generating confined conditions in this lower aquifer. This last aquifer receives groundwater from lateral and vertical flows from the carbonate hill of Mt. Petrino, which is not hydro-geologically connected with the Mt. Massico, and from the carbonate bedrock, both shallow and deeper groundwater flow from NE to SW, toward the sea [22,60] (Figure 2).

Groundwater below the tuff presents interesting chemical characteristics: temperatures and average electrical conductivities are above 20 °C and >20,000 µS/cm, respectively, while, in the neighboring areas, these parameters show lower values (i.e., conductivity around 1000 µS/cm). These anomalies can be explained by the rise of deep mineralized waters (rich in gas, above all CO₂ of inorganic origin [61]) along the faults of the carbonate bedrock. These deep waters mix with groundwater below the tuff and increase their total dissolved solids and temperature, features that are distributed along the groundwater flow direction [22].

Additional data about the mineralized zone is derived from drilling, in the frame of VIGOR project, of a geothermal exploration well, 300 m deep, about 150 m from the base of Mt. Petrino (Figure 2). The first 167 meters of the well stratigraphy (

Table 2. Geothermal stratigraphy and geophysical properties of the layers composing the analyzed domain.

Layer Number	Depth from Ground Level (m)	Lithology	Density (kg/m3)	Permeability (m2)	Porosity (-)	Thermal Conductivity (W/(m K))
1	0–1.20	Silty-clayey vegetal soil	1600	1.02 × 10−12	0.50	1.00
2	1.20–2.80	Gravel in silty-clayey matrix	1800	1.02 × 10−10	0.30	2.20
3	2.80–5.50	Clays and silts, with sandy and peaty intercalations	1450	1.02 × 10−14	0.60	1.50
4	5.50–15.0	Peat pits and clays with sandy and coarse interactions	1100	1.02 × 10−12	0.65	2.00
5	15.0–25.0	Grey tuff	2000	1.02 × 10−13	0.47	1.10
6	25.0–32.7	Sands alternated with clayey silts	1700	1.02 × 10−11	0.35	2.00
7	32.7–34.8	Limestone from fractured to compact	2000	1.00 × 10−13	0.15	2.80
8	34.8–51.4	Highly fractured and karst limestone	2200	1.00 × 10−12	0.15	2.80
9	51.5–101	Highly fractured and karst limestone	2100	1.00 × 10−11	0.30	3.00
10	101–123	Milonitic rock, sometimes cemented friction breccia	2500	1.00 × 10−10	0.15	2.80
11	123–157	Very fractured limestone with karstic voids filled with loose material	2100	1.00 × 10−11	0.15	2.80
12	157–167	Very fractured limestone, calcite and manganese concretions in the fractures, episodic empty karst. Strongly fractured and karstic	2200	1.00 × 10−12	0.15	2.80

) highlight that limestone starts under pyroclastic-alluvial deposits at about 30 m b.g.l., and that fracturing of crossed carbonate rocks reaches its maximum around 100–120 m b.g.l. (

Table 3. Temperatures and lithology of the geothermal exploration well at different depths.

Depth from Ground Level (m)	Temperature (°C)	Lithology
0–32	38.5–40.0	Pyroclastic–alluvial deposits
32–124	40.0	Limestone passing to Milonitic rock
124–160	40.0–36.0	Milonitic rock passing to very fractured limestone
160–240	36.0–35.0	Limestone passing to very fractured limestone
240–300	33.0–30.0	Very fractured limestone passing to fractured limestone

).

Table 2 reports the geothermal stratigraphy and the geophysical properties of the layers composing the analyzed domain. Density, thermal conductivity, and porosity were derived from literature, whereas the permeability was determined through “Lugeon tests” performed at different depths in the exploration well.

The experimental campaign has revealed the presence of artesian groundwater (under a maximum pressure of 1.70 bars) in the carbonate rocks, from which, in absence of confinement, about 23.5 L/s naturally flow [1]. After reaching the depth of 300 m, inside the well (and with a temporary casing) geophysical logs were performed (using a probe 2PFA-1000 / MATRIX converter) to measure temperature (Table 3) and verticality. During drilling, groundwater samples have been collected at different depths to analyze changes in chemical composition. The gases have been also sampled at two different sites, in correspondence of the geothermal exploration well and in its vicinity [60]. It was found that the chemical characteristics of groundwater along the depth are similar; the gases sampled are 98% CO₂, whose average concentration is 1380 mg/L [1,60].

Combining all data, a geothermal model that controls the mineralization of groundwater at the base of Mt. Petrino can be proposed [1,10,22,26,60,61]. According to this model, the geothermal reservoir corresponds to the carbonate bedrock of the Mondragone plain. Near Mt. Petrino, deep hot gases (mainly CO₂) rise along the faults of the bedrock, involving groundwater (sodium-bicarbonate type) typical of a reducing environment and of meteoric origin. From the carbonate bedrock, thermo-mineral groundwater spreads in the overlying pyroclastic-alluvial aquifer diluting gradually. Therefore, the geothermal process is closely linked to the faults, as shown by temperature increase, which occurs, as measured continuously along the well (Table 3), at the same levels of major fracturing degree. Therefore, it is very likely that these levels represent the preferential path along which the hot fluid rise occurs [60,61].

4. Numerical Simulation of Geothermal Exploitation

In order to investigate the sustainability of geothermal exploitation, a numerical model was developed to analyze the coupled multiphase thermo-hydraulic processes due to the extraction and reinjection of groundwater. The model was implemented within the Finite Element commercial software COMSOL Multiphysics®, which is a powerful tool to reproduce coupled or multiphysics phenomena.

4.1. Governing Equations

Geothermal energy exploitation involves the interaction of different physical phenomena occurring in the ground, such as fluid flow, heat transport, chemical transport, and mechanical deformation [32,34,62]. In this work, only heat and fluid flow in fully saturated porous media were analyzed, considering the conservation laws of mass, momentum, and energy in a porous medium, that represents the ground.

The conservation of mass for a fluid flowing through a porous medium is expressed as:

$$\frac{\partial }{\partial \mathrm{t}}\left(\rho \epsilon \right)+\nabla ·\left(\rho \mathbf{u}\right)=0,$$

(1)

where ρ is the fluid density, ε is the porosity, u is the seepage (Darcy) velocity vector. Therefore, the terms ρε and ρu represent the mass per unit volume within the porous matrix and the fluid mass flux, respectively. The flow velocity is low due to the low values of permeability and porosity of the domain. In this study, the porous medium flow was analyzed by using the well-known Darcy equation, which is the most common approach in geo-mechanic [63]. According to Darcy equation, the net flux across a face of porous surface, u, is linearly related to the pressure gradient, ∇p, as

$$\mathbf{u}=-\frac{\mathrm{K}}{\mu }\nabla \mathrm{p},$$

(2)

where K is the permeability of the porous medium and µ the geothermal fluid dynamic viscosity. Since the fluid flow is coupled to heat transfer, the dependence of fluid properties on temperature is taken into account.

The Darcy’s flow model is coupled with heat transport in a porous medium to study the temperature distribution under geothermal exploitation. In the heat transport equation, the local thermal equilibrium between the solid and fluid phases is considered. Therefore, the solid temperature, T_s, is equal to the fluid temperature, T_f. Moreover, heat transport between solid and fluid phases is considered to be negligible. The heat transport in the subsurface is described as [64,65]

$${\left({\rho \mathrm{c}}_{\mathrm{p}}\right)}_{\mathrm{eq}}\frac{\partial \mathrm{T}}{\partial \mathrm{t}}+\left({\rho \mathrm{c}}_{\mathrm{p}}\right)\mathbf{u}·\nabla \mathrm{T}=\nabla ·\left({\mathrm{k}}_{\mathrm{eq}}\nabla \mathrm{T}\right),$$

(3)

where c_p is the specific heat, k is the thermal conductivity. The term (ρc_p)_eq represents the volumetric heat capacity of the porous medium. The subscript eq indicates that these are equivalent properties since they are spatially averaged to account for the porous matrix according to

$${\left({\rho \mathrm{c}}_{\mathrm{p}}\right)}_{\mathrm{eq}}{=\theta }_{\mathrm{p}}{\rho }_{\mathrm{p}}{\mathrm{c}}_{\mathrm{p},\mathrm{p}}+\left({1-\theta }_{\mathrm{p}}\right){\rho \mathrm{c}}_{\mathrm{p}},$$

(4)

$${\mathrm{k}}_{\mathrm{eq}}{=\theta }_{\mathrm{p}}{\mathrm{k}}_{\mathrm{p}}+\left({1-\theta }_{\mathrm{p}}\right)\mathrm{k},$$

(5)

where θ_p is the porous matrix volume fraction and subscript p indicates quantities that are related to the solid porous matrix.

4.2. Boundary and Initial Conditions

The boundary conditions employed in the present model refer to the ground and the well domains. For logistical reasons, the position of the injection well is not the traditional one, which is downstream of the extraction well [66]. The injection well is in fact placed, along the direction of groundwater flow, parallel to the extraction well and not far behind this (Figure 2). The top and bottom surfaces of the ground domain are considered impermeable to mass (u = 0), such as the lateral surfaces parallel to the groundwater flow. The average hydraulic gradient, (i), in the groundwater layer was estimated, according to Corniello et al. [22], with a NE–SW direction of the flow, determined from the iso-piezometric map of the area (Figure 2).

A prescribed velocity u was assigned on the filtering lateral surfaces of the wells. This value was calculated considering the water flow rate of extraction and injection, $\stackrel{.}{V}$, measured during the experimental campaign. As initial condition for the groundwater flow equation, the hydraulic head, H, derived from experimental measurements, was imposed on the whole domain. The bottom ground surface was assumed to be adiabatic, whereas a convective heat flux is imposed on the top surface of the domain. The external temperature was assumed to be equal to 16.0 °C, which is the yearly average temperature of the area. As concerns the wells, a fixed temperature, T, equal to that of flow extraction and injection, is imposed on the lateral surfaces. The water temperature values measured along the exploration well [1,56] are illustrated in Figure 3. The soil was assumed to be in local thermal equilibrium with the groundwater, then the measured temperature profile (Table 3) was considered as the initial condition on the ground domain. All boundary and initial conditions are summarized in Figure 3 and

Table 4. Boundary and initial conditions.

Parameter	Symbol	Value	Unit
Hydraulic gradient	ix	0.5%	m/m
Extraction flow rate	${\stackrel{.}{V}}_{out}$	84.6	m3/h
Reinjection flow rate	${\stackrel{.}{V}}_{in}$	84.6	m3/h
Hydraulic head	H	13.0	m
Extraction Temperature	Tout	34.5	°C
Injection Temperature	Tin	32.0	°C
External Temperature	Text	16.0	°C
Convective heat transfer coefficient	h	25.0	W/(m2 K)

.

The experimental campaign revealed a sub-vertical fault zone, NW–SE oriented. The injection well is located at 135 m (Figure 2) from the production well, in the rearward direction. From the stratigraphical data of the exploration well, two different faults were found between 100–120 and 180–240 m b.g.l., according also to the results of the detailed geophysical survey. The deepest fault if projected towards the injection well may or may not meet this well at about 80 m depth. For this reason, two case studies were analyzed, one with the fault crossing the injection well and the other without. The fault zone presents the same characteristics of the layer identified by the number 10 in Table 2.

5. Results and Discussion

The numerical model described in the previous section was used to study the sustainable use of the low-temperature geothermal field of Mondragone. In order to evaluate the sustainability of this use, the temperature field in the ground was modeled over a period of time in which the source was continuously used at the full capacity of the plant. As mentioned above, the analysis was carried out with and without the fault zone crossing the injection well, to assess its influence on heat transfer and fluid flow. The simulations that refer to a single fault are indicated as Case 1 from now on, whereas those with the domain that includes two fault zones are referred to as Case 2. Figure 4 shows a schematic representation of the domain considered for the analysis of Case 1 where a single fault, indicated as F1, is present. In Case 2 the domain is the same, with the addition of a fault zone crossing the injection well, as shown in Figure 5. For the sake of clarity, only the data concerning the fault crossing the injection well, indicated as F2, are reported in Figure 5.

The overall domain is defined as a rectangular parallelepiped of 163 m in depth with a basis equal to 250 × 200 m. Below this depth, the effect of the geothermal system is considered negligible due to the lower degree of fracturing in limestone. The fault zones have both a thickness of 20.0 m and a slope of 10° in relation to the vertical axis. The production and reinjection wells cross the fault cataclastic zone starting from 80.0 to 120 and 68.0 to 110 m b.g.l., respectively. Both the wells have a radius of 24.4 cm and a filtering section of 10.0 m in height from the top of the cataclastic zones.

The water properties, ρ, µ, c_p, and k are considered temperature-dependent. For the ground domain, literature values of thermal conductivity, density, porosity, and permeability of the eight layers found during the well drilling were considered. The parameters used in the numerical simulation are reported in Table 4.

A 3D unstructured mesh, refined near all the boundaries to capture the larger gradients of both flow velocity and temperature, was used. Indeed, the grid needs to be finer in the region of greatest change in the skin zone close to the well and can be coarser further out in the reservoir [67]. For the sake of clarity, the mesh used for the case of a single fault is reported in Figure 6.

The details of the meshes used for the two cases are shown in

Table 5. Mesh details for Case 1 (single fault) and Case 2 (two faults).

Parameter	Case 1	Case 2
Tetrahedral elements	3,373,929	3,671,990
Triangular elements	168,494	182,727
Edge elements	6134	6813
Vertex elements	160	180
Minimum element quality	0.03311	0.03311
Average element quality	0.6594	0.6601
Element volume ratio	3.95 × 10−6	2.10 × 10−6

. A mesh sensitivity study was carried out to obtain grid-independent results from the calculations.

In order to assess the impact of the geothermal exploitation on the ground domain, the groundwater velocity field and the ground temperature field were determined through the proposed model for a period of time of 60 years.

The temperature fields at different times in the xz plan cutting the domain in correspondence of the two wells’ axes are reported in Figure 7 for the case of the single fault (Case 1) and in Figure 8 for the case of double (Case 2) fault zones. For the sake of clarity, the section cuts are also reported.

The effect of geothermal extraction is mainly contained in the area surrounding the production well. A moderate decrease in the ground temperature can be observed in the downstream region: the variation is similar in both Cases 1 and 2.

In the upstream section the temperature field remains almost unchanged in Case 1, whereas in Case 2 a temperature decrease can be observed in the region delimited by the two faults. The difference is due to the permeability and the thermal conductivity of the fault zone, which are larger than those of the surrounding material. However, in both cases after 10 years of operation the temperature field does not change any more (reaching steady state operation of the reservoir). The effect of injection is more pronounced than that of extraction in both cases. However, the thermal disturbance does not reach the superficial layers, since it is limited by the 5th layer which has lower values of permeability and thermal conductivity with respect to the other layers (Table 2). The injection in the fault slightly modifies the temperature distribution within the domain: the propagation of the thermal disturbance is lower in layers over the fault and higher towards the bottom of the domain. In both cases, the temperature field does not change any more after 20 years of operation, therefore the period of time needed to reach a steady condition is longer in the region surrounding the injection well than in the area of the production well. This is probably due to the higher difference in temperature between groundwater surrounding the injection well and injected water.

The temperature field, at different times, in the yz plan cutting the domain in correspondence of the axis of two wells is shown in Figure 9 and Figure 10 for Cases 1 and 2, respectively. For the sake of clarity, the section cuts are also reported. In these plots, it can be noticed that the temperature decrease caused by the injection well is localized in the region of the domain behind the production well and it is less pronounced in Case 1. In Case 1 the temperature propagation is more pronounced towards the upper layers, whereas in Case 2 it is larger in the bottom of the domain.

The temperature distribution and the streamlines indicating the velocity field, at different depths and times, reported in Figure 11, clearly show that the direction of the thermal disturbance is significantly influenced by the groundwater flow. A stable temperature distribution is reached after 25 years at the lowest depths, whereas only 5 years are needed at highest depths. The thermal disturbance in the layer close to the top surface of the domain is concentrated in the region surrounding the injection well.

The fault located in correspondence of the injection well, F2, does not significantly affect the temperature distribution, except for the upstream region of the production well, where the temperature is slightly lower than in Case 1. Conversely, the fault influences the velocity field, due to the permeability, which is larger than in the adjacent layers. This variation in the velocity field is mainly responsible for the extension of the region affected by the thermal disturbance in Case 2.

At depths below 75.0 m from ground level the temperature and velocity fields are constant over time, therefore only their variation with depth was taken into account. Analyzing the plots in Figure 12, which are all related to the time of 60 years, it can be observed that the thermal disturbance is more extended in Case 2. In addition, at 155 m b.g.l. the thermal disturbance can be neglected in Case 1, whereas in Case 2 the region downstream of the injection well is still affected by the geothermal exploitation.

The Darcy velocity field and magnitude at different depths for Cases 1 and 2 are reported in Figure 13 and Figure 14, respectively. The xy sections reported in these figures are related to different layers of the domain, in order to analyze the influence of their properties on the groundwater flow. In the first two layers of the domain, (refer to Figure 13a,b and Figure 14a,b), both the distribution and the magnitude of Darcy velocity do not change from Case 1 to Case 2, and result as constant in the entire xy plans considered. This is due to the low permeability of the third layer, which avoids the effect of the geothermal exploitation on the groundwater flow reaching the top of the domain. Indeed, the groundwater flow in the lower layers of the domain, which are characterized by higher values of permeability, varies in the regions where the production and injection wells are located (refer to Figure 13c,d and Figure 14c,d).

Considering the injection well, the velocity magnitude decreases upstream and increases downstream, whereas the opposite trend can be observed in the region surrounding the production well. In Case 2, a lower propagation of the fluid flow disturbance can be observed: this is due to higher permeability of the material in the correspondence of the fault, which facilitates the groundwater flow, which is more concentrated in the tectonic region. In addition, considering Case 2, the velocity magnitude varies differently in the zones of injection and production, probably due to the different vertical extension of the faults. The fault located in the injection area, F2, is more extended and its effect on the upper layers is larger.

In the sections that cut the wells (refer to Figure 13f,g and Figure 14f,g), the Darcy velocity field and magnitude are similar in Case 1 and Case 2. Obviously, the magnitude is significantly higher than at other depths due to the flow extraction and injection. A different trend can be observed in sections 10.0 m distant from the bottom of the wells in the vertical direction (refer to Figure 13e,h and Figure 14e,h), mainly in the areas surrounding the injection well. The fault affects the groundwater flow, whose magnitude increases in the area surrounding the injection well and decreases upstream of the production well. However, a vertical distance of 10.0 m is sufficient to observe a significant decrease in the Darcy velocity, which tends to the undisturbed value as the horizontal distance from the wells increases. At higher depths, the velocity magnitude continues to decrease (refer to Figure 13i,l and Figure 14i,l). The velocity magnitude is higher in the central region of the considered section, where the groundwater has the same direction of the undisturbed fluid flow in both cases, injection and extraction.

Finally, in Figure 15, the temperature over time is plotted at three different distances from the wells (10, 50, and 100 m) at a depth equal to half of the filtering sections of the wells, for the two cases analyzed.

Considering the production well (see Figure 15a,b) whose filtering section goes from −90 to −100 m, the geothermal exploitation causes a temperature decrease of around 2 °C during the considered period of time. The temperature becomes constant after 10 years in Case 1 and 20 years in Case 2. It is worth noticing that in Case 2, the temperature decrease at a distance of 100 m is slightly lower than in Case 1, suggesting a positive effect of the injection in the fault. A similar conclusion can be drawn analyzing the results related to the injection well (see Figure 15c,d). The variation in temperature with respect to the initial condition is lower than 1 °C, therefore significantly lower than that observed in the region of the production well. The period of time needed to reach a constant temperature is longer than for the production well: 20 years for the Case 1 and 25 years for the Case 2. The injection in the fault has a positive effect, reducing the temperature oscillations during the first years of exploitation in Case 2. The low decrease of the basin temperature is coherent with the results found in the literature related to the geothermal exploitation of low-temperature reservoirs [2].

So far, the numerical results show that the geothermal exploitation modifies the groundwater flow in correspondence with the filtering sections of the wells, whereas only a weak effect can be observed in the rest of the domain, in terms of speed and direction. In order to account for the possible presence of a fault crossing the injection well, two case studies were analyzed. When the injection occurs in the fault the variation of the velocity field is lower, thanks to the permeability in correspondence with the fault. The injection in the fault positively affects the temperature distribution as well. The thermal disturbance does not reach the top layers of the domain due to their hydraulic and thermal properties in both cases. It is worth noticing that in the case of injection in the fault the disturbance tends to be more concentrated in the same fault region, decreasing its effect on the other layers of the domain. In addition, in this case, also the temperature oscillation over time observed at the depth of the wells is lower than in the case without the injection in the fault.

The proposed model allows estimation of the effect of geothermal exploitation not only on the reservoir but also on the operation of the technological system that uses the geothermal source for thermal energy production. The results of the simulations indicate a temperature decrease of around 2.00 °C in the extraction temperature during the considered period of time. Such a decrease will reduce the thermal energy available from geothermal exploitation. In the specific case analyzed in this work, where a heat pump was considered to supply heat to two schools, a decrease in the COP of the heat pump is expected.

6. Conclusions

The numerical results on the geothermal exploitation here reported are of global interest for the following reasons: (a) the nature of aquifer represented by carbonate rocks with strong anisotropy, due to the presence of several fault systems, with confined groundwater; (b) the environmental effects of an open-loop well-doublet plant are assessed considering the geometry of the wells (imposed by logistical constraints) which is designed in a different way from the one usually adopted, where the injection well is placed downstream of the extraction well. In fact, even if thermal feedback is present at the production well due to the position of injection well, the sustainability of the use of the geothermal resource over time (60 years here tested) is asserted by the numerical simulation. In fact, the main results show that the temperature decrease is only 2.00 °C in correspondence with the production well filtering section, and even less at greater depths.

It is likely that the sustainability is linked, considering the hydrogeological scheme, to the characteristics of groundwater flow (direction and pressure) and to the position of the filtering section of the extraction well, that is placed near the main cataclastic zone of a fault along which the permeability is higher than the surrounding rocks and the upwelling of deep hot fluids is easier.

Even assuming the presence of another fault at the injection well, the numerical simulation still indicates the full sustainability of the geothermal plant. It is worth noticing that in the case of injection in the fault the disturbance tends to be more concentrated in the same fault region, decreasing its effect on the other layers of the domain. Finally, in all the cases considered, the thermal disturbance does not reach the top of the domain due to the hydraulic and thermal properties of the rocks closest to the ground level.

The definition of the hydrogeological and stratigraphic scheme of the geothermal area together with the numerical modelling are therefore decisive either for the final choice of the well plant localization or for the evaluation over time of the environment effects of the use of the geothermal resource in faulted carbonate rocks, which are common in Southern Italy and peri-Mediterranean regions.