The Castel Giorgio-Torre Alfina (CG-TA, central Italy) is a geothermal reservoir whose fluids are hosted in a carbonate formation at temperatures ranging between 120°C and 210°C. Data from deep wells suggest the existence of convective flow. We present the 3D numerical model of the CG-TA to simulate the undisturbed natural geothermal field and investigate the impacts of the exploitation process. The open source finite-element code OpenGeoSys is applied to solve the coupled systems of partial differential equations. The commercial software FEFLOW® is also used as additional numerical constraint. Calculated pressure and temperature have been calibrated against data from geothermal wells. The flow field displays multicellular convective patterns that cover the entire geothermal reservoir. The resulting thermal plumes protrude vertically over 3 km at Darcy velocity of about  m/s. The analysis of the exploitation process demonstrated the sustainability of a geothermal doublet for the development of a 5 MW pilot plant. The buoyant circulation within the geothermal system allows the reservoir to sustain a 50-year production at a flow rate of 1050 t/h. The distance of 2 km, between the production and reinjection wells, is sufficient to prevent any thermal breakthrough within the estimated operational lifetime. OGS and FELFOW results are qualitatively very similar with differences in peak velocities and temperatures. The case study provides valuable guidelines for future exploitation of the CG-TA deep geothermal reservoir.

1. Introduction

Since the 1970s, the increasing threat of a worldwide energy crisis has prompted many governments to reduce their dependence on traditional nonrenewable energy sources focusing on renewable ones (e.g., geothermal energy, hydroelectric, wind-energy, and several forms of solar energy, such as bioenergy, biofuel, photovoltaic, and solar-thermal energy [1]). Geothermal energy is expected to play an increasing role to meet future power demand. This is related to its enormous exploiting potential, which is capturing the attention of industries also due to technological advances in the exploration of promising geothermal fields [2, 3].

Most of the world famous geothermal energy sources exploited today are associated with volcanic and/or recent tectonically active areas. Important examples are the geothermal areas of Yellowstone [46], northern California [7], the Pannonic Basin [8], and the Rhine Graben [911]. Italy strongly contributes to the development of geothermal power generation. In 1904, the world first electrical power was produced from a geothermal energy source in the Larderello site (Tuscany, central Italy) [1, 1218].

Beyond the world famous Larderello system, fossil and active hydrothermal manifestations are distributed all along the Preappennine belt of central Italy, facing the Tyrrhenian coast. This area has undergone both lithospheric extension and upper mantle doming. Such processes have been active since the Miocene [19, 20] and are likely sustained by mass and heat fluxes from the upper mantle. This is suggested by the intense tectonic and volcanic activity associated with extremely high and variable surface heat-flux anomalies [21, 22]. All these processes document a predominant heat transfer mechanism by vertical mass flow, which accumulates large amount of geothermal resources at accessible depths in the upper crust. Two geothermal fields in the area, characterized by heat flow values of up to 1000 mW/m2 (Larderello system) and 600 mW/m2 (Mt. Amiata system) [11], are currently exploited for the production of electricity.

Recently, various projects were set up on a regional basis to investigate the geothermal potential of the Italian Tyrrhenian facing areas. Moreover, research and development of new exploitable geothermal fields have been encouraged by the approval of specific decrees of law (i.e., Legislative Decree of 11 February 2010, n. 22, modified by Legislative Decree of 3 March 2011, n. 28 and Article 28 of Decree of Law of 18 October 2012, n. 179).

Among the identified promising areas, the Castel Giorgio-Torre Alfina field (CG-TA, northern Latium, Figure 1) is an example of an early explored and so far not exploited medium-enthalpy geothermal system [2325]. Detailed hydrogeothermal data, available for the selected area since early 70s, show that the CG-TA is a potential geothermal reservoir with medium thermal characteristics (120°C–210°C) whose fluids (pressurized water and gas, mainly CO2) are hosted in a fractured carbonate formation [24, 2630]. Data from the deepest geothermal drilling in the area (Alfina015 well, max depth −4,826 m a.s.l., see Figure 1 for location) show a highly variable temperature gradient ranging between 0.15°C/10 m and 2.1°C/10 m [31]. Such a strong variation likely indicates the presence of highly convective flow within the reservoir rocks. This finding makes the CG-TA area suitable for future exploitation through a new generation 5 MWe geothermal pilot power plant. Following the guidelines of the above-mentioned Italian legislative decrees, this exploitation project is characterized by no gas emission to the atmosphere and total reinjection of the geothermal fluid in the same producing geological formation (i.e., geothermal well doublet system).

Planning of such challenging geothermal field exploitation projects requires an appropriate numerical modeling of the involved heat and fluid transfer processes. In the last 20–25 years, models have been set up for more than 150 geothermal fields worldwide [18, 22, 3243]. These numerical models allow us to define well’s system design, fracture paths, extraction rates, and temperature of injected and produced thermal waters, to interpret hydraulic tests or stimulation processes, and to predict reservoir behavior during geothermal power production. Therefore, they are mandatory to optimize the productive capacity and the thermal breakthrough occurrence [1, 44].

Mathematical modeling of a geothermal reservoir allows reconstructing both the deep natural fluid circulation and physical/chemical fluid characteristics. This can be of interest at geothermal sites where high temperatures and strong corrosion, caused by very acidic involved fluids, occur. In some cases, the fluids may react chemically with the hosting rocks, precipitating minerals that diminish reservoir permeability by pores and fractures obstruction [45]. These phenomena create spatially variable patterns of mineralization and permeability, thus affecting the exploitation of the reservoir [45, 46].

Numerical modeling of exploited geothermal systems should include (i) a solid conceptual model of the reservoir geology and structure, (ii) the location and geometry of wells and possible fractures systems, and (iii) the parameterization of hydraulic, thermal, mechanical, and chemical (HTMC) properties of the reservoir and of the involved fluids [47].

The aim of the present study is to build the 3D numerical model of the deep, medium-enthalpy CG-TA reservoir to reproduce the highly convective undisturbed present-day natural state of the reservoir. These results, validated against the pressures and temperatures measured in geothermal wells, are afterward used to investigate the feasibility of a geothermal power production configuration (i.e., injection and production wells). The analysis is performed on a hypothetical 50-year operational life cycle adopting a well doublet system at a 1,050 t/h flow rate [48]. The finite-element open source code OpenGeoSys [49] is used to build the hydrothermal (HT) model. As additional numerical constraint, the results are compared against those obtained with the commercial finite-element code FEFLOW [50].

First, the hydrogeothermal data derived from geophysical investigations and from geothermal wells are described and used to build a conceptual and numerical model of the CG-TA reservoir. Then, the numerical approach based on the OpenGeoSys software is given. Results are obtained both at short-term (i.e., operational) and long-term (i.e., full reservoir recovery) time scale. Besides providing valuable guidelines for future exploitation of the CG-TA deep geothermal reservoir, this study highlights the importance of field data constraints for the interpretation of numerical results of fluid processes in reservoir-scale systems.

2. Reservoir Characterization

2.1. Regional Geological Setting

The occurrence of medium- and high-enthalpy geothermal fields in central Italy is localized along the Tyrrhenian margin of the Apennines (Figure 1). The complex geologic and tectonic settings of this area have been studied by several authors ([5154] and references therein).

The present-day structural setting of the Tyrrhenian coast facing regions represents a heritage of compressive and extensional geodynamic processes that began in the Oligocene (i.e., 30 Ma BP) with the Alpine-Apennine orogenesis [18, 55]. The compressive phase resulted in the formation of fold-and-thrust-belts and associated piggy-back basins with NNE-SSW oriented trend [5658]. Then, the subsequent extensional phase, due to the Tyrrhenian back-arc extension, resulted in the formation of NW–SE tectonic basins and in the crustal thinning, with consequent upwelling of magma bodies and increased heat flow [18, 21, 22, 59, 60].

Due to the interplay of all these phenomena, the geologic and structural settings of the area are quite complex and involve many different lithostratigraphic units. The main and most widespread complexes, from the shallower to the deeper ones [18, 25, 61, 62] (Figure 2), are, namely,(i)volcanic complex: volcanic products including tuff, lavas, and pyroclastic rocks, characterized by variable thickness, with a maximum of 200 meters,(ii)neoautochthonous complex: clays, with limited sand content, conglomerates, and detrital limestones in a discontinuous layer 50 to 160 meters thick,(iii)Ligurian/sub-Ligurian complex: Jurassic-Eocene clayey-marly units in flysch facies, sandstones, marly limestones, and ophiolites. They are characterized by a highly variable thickness ranging from 500 to 1,800 meters (RAI01 well, see Figure 1) [61],(iv)Tuscan and Umbria Nappe complex: Triassic-Lower Miocene arenaceous and clayey-marly formations, calcareous-siliceous rocks, dolostone, and anhydrites. The upper portion of this formation is characterized mainly by marly limestone and shales and is referred to as the “Scaglia formation.” The Tuscan and Umbria Nappe carbonatic formation reaches a thickness of about 3,700 meters (Alfina015 well, see Figure 1) [31].

2.2. The Castel Giorgio-Torre Alfina Geothermal Field

The CG-TA geothermal field (Figure 1) is located to the north of the Vulsini caldera [24], at the boundary between the Tuscany, Umbria, and Latium regions (central Italy). The Torre Alfina reservoir was extensively explored between the 1970s and the 1990s. We refer to the works of Cataldi and Rendina [23] and of Buonasorte et al. [24, 31] for the detailed description of the geothermal explorations carried out in the area. These investigations culminated with the drilling of eight geothermal wells, with depths ranging from 563 to 2,710 m, and more recently, with the drilling of a very deep geothermal well (Alfina015, Figure 2) reaching the depth of 4,826 m.

The integration between stratigraphic borehole logs and geophysical [24, 31] and seismic [27] data identified the CG-TA geothermal reservoir as hosted in a structural high (i.e., horst structure highlighted in the correlation section of Figure 2) of fractured Mesozoic limestones, belonging to the Tuscan and Umbria Nappe complex, and marked by positive geothermal and magnetic anomalies [63]. Structural investigations performed in the area by Buonasorte et al. [56] and Piscopo et al. [64] provide a detailed description of the N-S striking postorogenic extensional faults bounding this horst structure and an analysis of the geometry, orientation, and kinematics of all the other tectonic features occurring in the Torre Alfina geothermal system.

The first geothermal drilling campaign performed in the CG-TA field (1971-1972) was aimed to reach and cross the argillaceous and shaly terrains of the Ligurian and Sub-Ligurian complex. These investigations allowed not only a detailed stratigraphic reconstruction, but also the definition of the basic characteristics of the geothermal reservoir fluids (e.g., pressurized hot water with average temperature of 140°C) and the detection of a gas cap made by 2% of dissolved CO2. This 100-meter thick cap, recognized only in the central part of the field, was extensively exploited until few years ago for CO2 storage by the well Alfina013 (Figure 1).

The target of the more recent campaign (1987-1988) was the deeper and hotter geothermal reservoir, hosted in the metamorphic rocks lying underneath the calcareous formations. Though the exploration did not reach the metamorphic basement, it demonstrated the presence of a single very thick carbonatic reservoir (>3,700 m thick), within which a highly variable temperature gradient of 0.15°C/10 m–0.45°C/10 m was recorded [31]. These exploration wells resulted in multiple pressure and temperature vertical profiles within the geothermal field, three of which are illustrated in Figure 3. The available data stands in different depth ranges: Alfina002 well, the shallower one, with measured temperature data reaching −500 m a.s.l.; the second well (RAI01) reached −2,000 m a.s.l., while the last and most recently drilled Alfina015 well provided a full temperature profile up to a depth of −4000 m a.s.l. The shallower Alfina002 and RAI01 wells, reaching only the top of the reservoir units, registered a linearly increasing temperature with a high geothermal gradient in the range of 1.7–2.1°C/10 m [24]. This suggests a mainly conductive heat transfer mechanism associated with the cap rock impermeable units. A similar trend was observed in the shallower portion (up to ca. 1000 m depth) of the deep Alfina015 well (Figure 3). At about 1000 m a knick-point and a thermal inversion are observed along the profile. The geothermal gradient below this depth ranges between 0.15°C/10 m and 0.45°C/10 m. Such a strong variation, coupled with measured top and bottom nearly constant temperatures of the reservoir fluids (i.e., 140°C and 207°C at 1,050 m and 4,000 m depth, resp. [31]), point toward an intense, large-scale convective flow confined in the area of the buried structural high.

In summary, the CG-TA area is an example of a promising, early explored and yet to be developed geothermal filed. Despite the extremely favorable conditions for exploitation [23], its industrial development was not promoted till 2011. A new geothermal research permit was requested for the Torre Alfina area, aimed to the development of 2 new generation 5 MWe pilot doublet plants, with reduced gas emission [30, 48].

2.3. Conceptual Model

The spatial extent of the model is fundamental for a reliable simulation of the complex processes involved in a geothermal reservoir. An overly restricted scenario hampers a complete representation of the circulation into the field, whereas a very large one results in a more uncertain geological reconstruction and excessive computational loading. The model area, covering about 293 km2 (Figure 4), is located north of the Vulsini calderas and it is bordered by the Meso-Cenozoic ridge of the Mount Cetona, to the north and by the Bolsena caldera structure, to the south (Figure 1). The extent is large enough for the imposed boundary conditions not to interfere with the phenomena occurring inside the geothermal field. This is guaranteed by a horizontal distance between lateral model boundaries and the geothermal field of 7.5 km in the E-W direction and of 2.5 km in the N-S direction (Figure 4). Due to the large extent of the geothermal reservoir and its intrinsic geological complexity, a complete review of the existing data and literature was required [23, 25, 30, 31, 57, 61, 63, 65].

The geological model was based on deep geological cross sections [24, 61, 65] and contour line maps of the contact surfaces between geological formations [24, 25]. Major attention was devoted in representing changes inside and outside the geothermal reservoir. The base of the model was located at −4,500 m a.s.l., within fractured limestone reservoir units. The upper limit was defined by a rather flat topography derived from a 20 × 20-meter DEM derived topography. This resulted in a maximum model thickness of about 5 kilometers (from +600 m a.s.l. to −4,500 m a.s.l.; see Figure 4).

The reservoir units are composed, from bottom to top, of evaporites, limestones, marls, and radiolarites (Tuscan and Umbria series s.l. [65]). Such reservoir units are buried by the sealing units, and crop out at San Casciano dei Bagni village. The sealing units consist of an allochthonous flysch-type sequence composed of arenaceous turbidites intercalated with layers of shales, marls, and limestones, overlaid by an ophiolitic sequence (siliceous shales and sandstones including blocks of gabbro and serpentinite) (Ligurian units s.l. [66]).

As previously mentioned, the area has undergone a strong postorogenic deformation phase, resulting in strike-slip and subordinate normal fault systems (with associated fracture network) cutting and dislocating the internal architecture of the reservoir [60]. No anomalous soil CO2 flux was recorded by the detailed investigations performed by Carapezza et al. [30]. This indicates the effectiveness of the impervious behavior of both the sealing units, which are continuous all over the reservoir area with a thickness of no less than 400 meters and the fault system connecting the geothermal reservoir with the surface.

In summary, the conceptual model consists of seven hydrogeological units (Figure 4), of which the upper three form the sealing cap and the remaining comprise the reservoir. The volcanic complex (1) is the youngest one and it outcrops only in the southern part of the model domain. This formation tends to thin towards the north where it is in contact with the other sealing units, represented by the neoautochthonous complex (2) and the Ligurian/sub-Ligurian complex (3). The shallower geothermal reservoir unit is referred as the Scaglia formation (4), a tiny layer mainly consisting of argillites. Below this, the fractured limestone rocks of the Tuscan limestone formation (5) and the deeper Umbria limestone formation (6) are emplaced. The Scaglia, Tuscan, and Umbria complex units (numbers 4, 5, and 6) were additionally subdivided between formations stacked into the proper geothermal reservoir (i.e., the real portion affected by convection phenomena, with an extent of ca. 73 km2) and those falling outside the producing area. The model includes also a NE-SW trending subvertical fault (7), with a surface trace of about one kilometer, a vertical extent of 1.5 km, and impervious behavior.

3. Numerical Modeling

3.1. Modeling Approach

Based on the conceptual model, a refined reservoir-scale three-dimensional thermohydraulic (TH) model was built to investigate the different processes involved in the CG-TA geothermal reservoir.

The open source finite-element simulator OpenGeoSys (OGS) [49] was used to solve the differential equations governing density-driven flows. The mathematical and numerical formulation of the problem and the strongly coupled system of equations can be found in Kolditz et al. [49]. OGS fully implements several equations of state (EOS) in order to reproduce temperature and pressure dependent fluid density and viscosity. Here we used the polynomial fittings introduced by Magri et al. [67] that are valid for a wide range of temperatures (°C) and pressures ( MPa).

The model surface was discretized into 17,768 triangular finite elements, satisfying Delaunay’s criterion, by using the GMS software [68]. Mesh refinement was applied to ensure simulation robustness: elements size decreases gradually from 500 meters, at model lateral boundaries, to 10 meters close to the fault zone and around the geothermal wells (Figure 5). We verified that a finer mesh did not affect the calculated patterns.

The 2D surface grid was extruded vertically using a fully unstructured tetrahedral 3D mesh. The total volume of the model was discretized with 35 layers ranging in thickness from 250 meters, at the model bottom, to a minimum of 10 meters, near the topographic surface. In total, the 3D mesh consists of 1,720,774 tetrahedral elements (Figure 5) that preserves all outcropping and internal pinching of the geologic formations.

The two modeling challenges are (i) recreating the present-day, highly convective, unexploited, natural state of the CG-TA geothermal system and (ii) performing the predictive analysis of the industrial exploitation process of the field. Two scenarios are therefore presented [25]. (1) The first one, referred henceforth to as “natural state simulation,” reproduces the thermohydraulic dynamic conditions of the geothermal reservoir, without extraction or injection of fluid. Pressure and temperature values measured in the three geothermal wells drilled in the area (Figure 3) were used to constrain the numerical results. (2) Once a qualitatively satisfactory match between calculated and observed patterns in these three geothermal wells was obtained, the calculated temperature and pressure fields were used to initialize the second simulation step. The latter includes the operating conditions based on a reasonable configuration of injection and production wells. This scenario, referred to as “exploitation process simulation,” also assesses the impacts of the exploitation process on the long-term (i.e., up to 10,000 years) natural geothermal flow of the reservoir after the production stage.

The same modeling framework (i.e., boundary conditions, initial conditions, equations of state, and spatial and temporal discretization) is applied to the finite-element commercial software FEFLOW.

3.2. Boundary Conditions

Temperature and pressure boundary conditions are summarized in Figure 5. In both scenarios, temperature and pressure distributions at the top were assumed to be time invariant. A fixed value of 15°C (i.e., Dirichlet type), corresponding to the average annual temperature of the area, and an atmospheric pressure value of 1 bar (i.e., Dirichlet type) were set. The implicit assumption is that the groundwater table and the ground surface coincide [23, 24]. Outside the reservoir area, temperature and pressure at the bottom boundary nodes were fixed too (i.e., Dirichlet type). The chosen values were calculated according to the average geothermal and pressure gradients of 0.3°C/10 m and 1 bar/10 m, respectively (Figure 5). On the other hand, given the anomalous geothermal gradient (1.7–2.1°C/10 m [24]) in the area of the buried structural high, an incoming heat-flux of 0.256 W/m2 (i.e., Neumann type) was applied at the nodes on bottom boundary below the reservoir area (Figure 5).

A no-mass flow condition was imposed over all the lateral boundaries (i.e., adiabatic and impermeable boundaries). As said above, the large distance between the grid boundaries and the reservoir area guarantees that applied boundary conditions do not affect the field behavior.

The “natural state simulation” was performed to determine the present-day reservoir condition, without any fluid extraction/injection scenarios. To let the system reach the present-day anomalous temperature field, the simulation covers a period of 1 million years. To verify the “natural state simulation,” the spatial distribution of the simulated temperature was compared with the measured thermometric vertical profiles in correspondence to 3 geothermal wells (Alfina002, Alfina015, and RAI01; see Figure 1).

To simulate field production and to predict the future system evolution, pressure and temperature boundary conditions remained those applied for the “natural state simulation.” A reasonable configuration of 5 production and 4 injection wells, separated horizontally by a distance of ca. 2 km [24], was inserted in the “exploitation process simulation” model (see Figures 1 and 5). A hypothetical 50-year production and injection time span, with a flow rate of 1050 t/h, was chosen following Buonasorte et al. [24], Marini et al. [69], and Colucci and Guandalini [25]. Starting from this production scenario, a flow rate of 210 t/h for each production well was applied. At each injection well, a constant injection temperature (i.e., Dirichlet type boundary condition) of 80°C and a 262.5 t/h injection rate were applied [48]. These boundary conditions, distributed over the nodes of the active length of the production/injection wells (ca. 300 meters discretized with 12 nodes), were set as time-dependent. At the end of the 50-year simulation run, the wells boundary conditions were removed and the simulation ran for an additional 10,000 years to investigate the recovery time and to test the technical sustainability of geothermal power production.

3.3. Initial Conditions

In preparation for the dynamic reservoir simulation, initial reservoir conditions have to be determined. These initial conditions include both the geothermal gradient and the fluid pressure gradient (i.e., when advection/convection is not involved). For this purpose, this initialization phase of the natural state was modeled as a steady state condition without the incoming heat-flux at the bottom of the reservoir. The temperature and pressure boundary conditions and model internal partitioning were set as described above. In this steady state initialization, temperature and pressure effects on fluid density and viscosity were neglected.

The values of the petrophysical parameters of the involved lithostratigraphic units (Table 1) were derived from available data for the area [2325, 30, 31, 57, 6163, 65, 70]. Default values for thermal conductivity of water (0.65 W/m/K) and heat capacity of water (4.2 MJ/m3/K) were used. This initialization resulted in a temperature field with values ranging from 15°C to 160°C, at the ground surface and the bottom boundary, respectively (Figure 6(a)). Fluid pressure ranges from 1 bar, at the ground surface, to 491 bar at the model bottom.

4. Natural State Simulation

4.1. Model Definition

The natural state simulation aimed to define the present-day, unexploited thermofluid dynamic conditions inside the geothermal reservoir including the advective/convective fluid motion. Simulation started by applying initial conditions defined as above. As common practice, the natural state simulations of geothermal fields require a long simulation time so as to attain pressure and temperature stabilization in the reservoir [18, 66, 71, 72]. Therefore, a 1 Ma simulation time has been chosen, neglecting effects of past climate change or transient effects in the rocks, and representing a generic geologic period. The performed transient simulation adopted a maximum time step size of 500 years. This time step coincides with the one used to update the fluid density and viscosity values as a function of calculated pressure and temperature. Applied boundary conditions remained the same presented above.

Regarding the assignment of required petrophysical parameters, the reservoir units were distinguished between those falling inside and outside the proper geothermal reservoir (see Figure 4). The latter were initialized with a value of permeability equal to  m2, while the remaining reservoir units preserved their typical fractured limestone permeability values, derived from the literature works and ranging from 10–14 to 10–15 m2 (Table 1). A very low permeability value of 10–18 m2 adopted for the overlaying sealing units allowed modeling their impervious behavior. A compressibility value of 10-10 Pa−1, slightly lower compared the one used for the reservoir units (i.e.,  Pa−1; Table 1), was assigned to these formations. This reduced compressibility allows maintaining the fluid pressure as simulated in the previous stationary system initialization phase, avoiding a fluid pressure rise due to temperature increase. The complete set of applied hydraulic and thermal parameters is given in Table 1.

4.2. Results of the Natural State Simulation

Results of the 3D convective flows are shown in Figure 6 for three different simulation times (i.e., 0, 20,000, and 125,000 years). From the initial conductive temperature field, equal to the average geothermal gradient of 0.3°C/10 m (Figure 6(a)), a very efficient convective circulation develops only into the geothermal reservoir units ((4b), (5b), and (6b), as listed in Table 1). This resulted in a gradual increase of temperature values in this area, while outside the producing units the pressure and temperature fields showed a full correspondence to those obtained at the end of the previous stationary system initialization. Fluid circulates in form of rolls and exhibits multicellular convective patterns, which start oscillating after ca. 20,000 years of simulation time (Figure 6(b)). This implies sharper inherent gradients and continuous creation and disappearance of convective plumes patterns. Within the producing area, three elongated convective cells stretched over the entire geothermal reservoir (Figure 6(b)). The observed cellular motion consists of multiple central up-flows, with a fluid velocity in the range of 2–4 10-8 m/s, and associated lateral down-flows. The strong convective behavior allows cold infiltrating groundwater to reach basement depths where it gets heated before starting its upward migration to the top of the geothermal reservoir. Comparing the results of the natural state simulation along 1D profiles with real surveyed thermal profiles, it is possible to identify the time instant for which the model fits the real reservoir conditions. The identification of the best-fitting simulated temperature profiles was performed through an iterative manual process by comparing computed 1D profiles against temperature profiles at 3 geothermal wells (Alfina002, RAI01, and Alfina015; see Figures 1 and 3). The attained best fitting occurs at the 125,000-year simulation time (Figure 6(c)) for all the three geothermal wells. At that time, the pattern of the three elongated convective cells is highlighted by a sharp difference in temperature between raising and sinking fluids. Velocity of the modeled convective cells rises to ca.  m/s, while maximum fluid temperature reaches 263°C.

The sigmoidal shape of the temperature profiles suggests the occurrence of the highly convective flow (Figure 7). In fact, all the three simulated profiles exhibit a clear thermal inversion as soon as the reservoir depth is reached (ca. −500 m a.s.l. for the Alfina002 well, ca. −1000 m a.s.l. for the Alfina015 well, and ca. −2000 m a.s.l. for the RAI01 well). In the upper 2 km of the temperature profiles (°C, conductive regime), the difference between simulated and measured values is at maximum 10°C (see Figure 7). A comparison between simulated and measured deep reservoir temperature values (°C, convective regime) was possible only for the Alfina015 well, representative of almost the entire thickness of the model. The computed profile, in correspondence to the Alfina015 well position, shows a well-developed trend with an almost constant temperature down to about −3,500 m a.s.l. This can only be associated with the presence of a convective cell. A good fitting of the Alfina015 profile temperature was obtained by slightly shifting the sampling profile location of a few meters so that it hits upward buoyant flow. FEFLOW and OGS models exhibit similar temperature-depth profiles in all wells, but at different simulation time (125,000 versus 230,000 years, see Discussion). Differences in temperature values are observed at maximum depth, where convection is dominant and controls the thermal evolution of the system (Figure 7).

The computed best-fitting natural state temperature field formed the initial condition for the following dynamic reservoir simulations of the effects induced by the production process.

5. Exploitation Process Simulation

5.1. Model Definition

Once a satisfactory match for the natural state is accomplished, a realistic scenario was set up for the future exploitation of the CG-TA geothermal field through a 5 MWe pilot doublet power plant.

To achieve this, the chosen configuration (see Figures 1 and 5 [24]) of 5 production wells (CG1, CG1A, CG2, CG3, and CG3A) and 4 injection wells (CG14, CG14A, CG14B, and CG14C) was inserted in the model. The production wells extract the geothermal fluids from the uppermost portion (from −300 to −700 m a.s.l.) of the reservoir units. The extraction depth range depends on the well position relatively to the top of the producing area. The injection of the 80°C fluid, at the above described rate, was designed at a depth ranging between −1,350 and −1,550 m a.s.l. For both production and injection sites, the well active length was fixed at 300 meters.

Pressure and temperature boundary conditions were those described for the exploitation process simulation (see boundary conditions, Figure 5). Initial conditions, mimicking the present-day distribution of both temperature and pressure field, were imported from the calculated best-fitting time step (i.e., 125,000 years of simulation time) of the previously performed natural state simulation. The exploitation time lasted 50 years and total simulation time of 10,000 years, so as to evaluate the field long-term effects induced by production process.

5.2. Field Exploitation

The evolution of well pressure over time, computed at a node with depth close to the well bottom, showed that the maximum differences, relatively to the initial pressure field, were reached at the end of the production time (i.e., after 50 years of simulation). In more detail, as shown in Figure 8(a) for wells CG2 and CG3, the production wells realized a depressurization in the 15–17 bar range, at the end of the first year of simulation, and then stabilized to an averaged value of 19 bar at the end of the production time. This pressure variation corresponds to approximately 12–14% of the initial pressure values (from 120 to 150 bar depending on the considered well) in the production wells. At the end of the exploitation of the geothermal field (i.e., after 50 years), no further fluid extraction occurred and the production wells exhibited a fast recovery. The monitored pressures raised back to the initial values in less than 100 years for the whole production site (see Figure 8(c) for the wells CG2 and CG3).

On the other hand, at the injection site, injected water resulted in strong overpressures rising quickly in the first years (i.e., around 14–16 bar) to stabilize to an average value of 20 bar after 50 years of simulation (see Figure 8(b) for wells CG14A and CG14B). These overpressures correspond to approximately 7–10% of the initial pressure field (from 225 to 240 bar, for all the wells) recorded in the injection wells. At the end of the production time, in the same way as for the production wells, pressure values recovered quickly to the initial undisturbed ones (see Figure 8(d)). Therefore, comparing model pressure distribution at the beginning (natural state) and at the end of the simulation (i.e., after 10,000 years), no significant variations could be observed.

The evolution of temperature over time for both production and injection wells is plotted in Figure 9. During system exploitation, the recorded temperature at the production wells exhibited a progressive increase over time (see Figure 9(a) for wells CG2 and CG1A). The difference between onset and end of production (i.e., after 50 years) temperatures varies from 2.5°C for CG2 well to a maximum value of 9.5°C for CG1 well. This increase in temperature resulted from the direct extraction of fluids from within a very strong convective system in the inner portion of the reservoir. Hence, from a thermal point of view, this analysis showed no interference effects between injection and production sites.

At the end of the phase of exploitation (Figure 9(c)), the recorded temperatures at production wells exhibited two different behaviors depending on the position of the well relatively to the generated convective cells. For example, in the CG1A well, the simulated temperature slowly decreased after the first 50 years of simulation time and recovered the initial undisturbed values in about 1,000 years. By contrast, the recorded temperature in the CG2 well firstly decreased to the initial value and then followed a gently increasing trend (ca. 2°C at the end of 10,000 years simulation, see Figure 9(c)). To investigate further this behavior, simulation time was extended to 30,000 years (Figure 9(e)). It turned out that, for wells located close to the convective cell (e.g., CG2 well), the recorded temperature exhibited strong convective oscillations, starting after 5,000 years from the end of production, and preventing the stabilization to the initial temperature values (see Figure 9(e)). This behavior was related to the evolution of the convection regime and to the pattern of multiple positive thermal anomalies.

Fluid temperature in the injection area slowly decreased over time, reaching the injected value of 80°C at the end of the production time (see Figure 9(b) for wells CG14 and CG14C). At the end of exploitation, after 50 years of simulation time, temperatures in the surrounding of the injected wells recovered the initial values after ca. 2,000–3,000 years (see Figure 9(c)).

To evaluate the thermal response of the CG-TA reservoir to the production process, the influence area of the “cold-water” front was investigated (Figure 10). At the end of the production time, the 80°C isosurface around the four injection wells covered a subspherical volume with ca. 1 km in diameter (Figure 10(b)). Therefore, the tested horizontal distance of about 2 km, between the production and injection sites, fully excluded the hypothesis of a thermal breakthrough.

6. Discussion

Within the present work, we set up an accurate hydrothermal model to recreate the highly convective behavior of the CG-TA reservoir and then simulate the exploitation of this undeveloped geothermal field. A general procedure for model calibration was applied [25, 35, 7375], consisting of a natural state modeling followed by an exploitation process simulation.

The above described natural state simulation resulted in a good match between simulated (via the OGS code) and measured temperature profiles after ca. 125,000 years of simulation time (see Figures 6(c) and 7). To analyze the temperature field at the best-fitting simulation time (125,000 years, Figure 11), the OGS model was sliced with a vertical plane along the E-W axis and passing through the Alfina015 well (A-A′ section in Figures 11(a), 11(b), and 11(c)). Two convective plumes were recognized within the reservoir (Figure 11(c)). Different vertical temperature profiles were extracted along this cross-section plane at different relative positions with respect to the convective plumes (i.e., interplume or in-between the two interacting plumes; axial or along the major plume axis; intersecting or crossing the upper plume head; outside or in a portion not strongly affected by a convecting plume, Figure 11(d)). An increase in the conductive temperature field of the sealing units (i.e., profile portion above the cap rock/reservoir contact elevation, ranging between −500 and −1,000 m a.s.l. depending on well position) is observed moving toward the axis of the plume (Figures 11(c) and 11(d)). Within the reservoir, a thermal inversion characterizes the temperature profiles which cut the main plume laterally (i.e., interplume, intersecting, and Alfina015 profiles). The sigmoidal shape of the temperature vertical profiles is similar to the one observed in many other convection-dominated geothermal systems, for which comparable analyses were performed [18, 22, 43, 47, 76, 77].

Multiple horizontal temperature profiles along A-A′ (Figure 11(e)) cross section were extracted at different depths (i.e., −4000, −3000, −2000, −1000, −500, 0, and 300 m a.s.l.). The deepest profiles clearly showed two positive thermal anomalies, with values reaching more than 200°C in correspondence to the plumes axis. The shallower profiles show a progressive merging of the two plumes.

Once the unexploited present-day temperature and pressure fields are determined, the computed natural state was used as initial condition to simulate field production and the future system evolution. The highly convective behavior of the system was suggested by the temperature graphs of the production wells in Figure 9(a). From these results, we conclude that the thermal breakthrough was prevented, as testified by the progressive increase in the recorded production site temperatures during the exploitation simulation. Moreover, the exploitation process induced only very small long-term changes with respect to the natural state of the geothermal system. In fact, at the end of the production time (i.e., after 50 years of simulation time), temperature in the production wells located close to the convective cell (e.g., CG2 well) exhibits strong convective oscillations, following the unexploited behavior modeled in the natural state simulation. Darcy velocity of such convective cells stands in the range of 7.5–8.5 10-8 m/s, therefore close to the preexploitation one.

The performance of the OGS code at modeling the convective flow within the geothermal system has been tested against the FEFLOW code. The two codes implemented the same equations of state. OGS and FEFLOW results show that the calculated patterns were qualitatively similar (e.g., multicellular convective fluid motion and velocities of convective cells), while differences in the calculated values existed (e.g., best-fitting time step in the natural state simulation, absolute pressure, and temperature values during exploitation).

Results of the FEFLOW and OGS simulation are plotted in Figure 7. The iterative manual identification of the best-fitting time step resulted in a good matching between FEFLOW simulated and real thermometric data around 230,000 years of simulation time. Even if a quite large time gap characterized the reservoir present-day situation modeled with the two software programs (i.e., best-fitting at 125,000 years of simulation time for OGS and 230,000 for FEFLOW), the simulated vertical profiles perfectly overlapped for the entire depth of Alfina002 and Alfina015 wells. As for the RAI01 well, a good match between the two tested software programs is observed in the shallower portion of the thermal logs (i.e., cover and impermeable units and conductive pattern). As soon as the reservoir depth is reached (i.e., −2000 m a.s.l.), convection is dominant. As a result, the simulated patterns can highly oscillate leading to larger temperature differences at selected simulation time steps (ca. 50°C at model bottom). At the best-fitting simulation time, in the shallower portion of the vertical temperature profiles (°C), the difference between real measured data and FEFLOW simulated values stands in the range of 5°C (see Figure 7). This difference increases in the deeper portion of the temperature logs ( > 150°C) due to the highly convective flow, as previously explained. After 230,000 years of simulation time, FEFLOW convective cells exhibited a maximum velocity of ca.  m/s and temperature values reaching 280°C.

Starting from the present-day unexploited temperature and pressure fields, the same production scenario was simulated with FEFLOW. The results of pressure and temperature versus simulation time for both production and injection wells are plotted in Figure 12. FEFLOW returned a trend very close to the one by OGS, for both pressure and temperature time evolution. It is worth pointing out that the gap in the pressure values (~15 bar, Figures 12(b) and 12(d)) is due to the fact that the two codes started from slightly different initial pressure fields, and thus the recovery process stabilizes to these initial undisturbed values. Furthermore, the applied initial pressure and temperature fields were defined on the natural state condition identified only by the available thermal logs. Therefore, model constrains were only applied to temperature while missing any present-day pressure data. As the two software programs started from the same initial temperature values identified in the natural state simulation, the time evolution during the field exploitation process is perfectly overlapped (~1°C gap, see Figure 12(e)). Moreover, at the end of field production, FEFLOW exhibits the same convective oscillations in the productive wells as already observed in the OGS results (Figure 12(g)).

Finally, the areal extent of the 80°C “cold-water” front was evaluated in FEFLOW as in OGS. The 80°C isosurface propagated away from the injection wells and reached its maximum extent at the end of the production time (i.e., 50 years). Again, a slightly irregular spherical shape, ca. 1 km in diameter, was observed. This confirms that the tested exploitation scenario prevents the thermal breakthrough in the same way as shown by OGS simulation.

7. Conclusions

The objectives of this study are to model the origin of the thermal anomaly observed in the CG-TA medium-enthalpy geothermal field, to investigate the feasibility of geothermal exploitation, and to test capabilities of different codes at modeling highly buoyant flows. A fit-for-purpose 3D numerical model of the CG-TA geothermal system was built using the open source OpenGeoSys (OGS) code and the commercial FEFLOW code. Following a general procedure for geothermal numerical models calibration, the present-day, highly convective, unexploited (natural) state model preceded the simulation of field production process. Starting from a steady state initialization of the reservoir, a satisfactory natural state modeling was achieved with limited differences between measured and computed temperatures. At higher depths, as convection is dominant, strong measured/calculated temperature discrepancies can be observed. The multicellular highly convective behavior (Darcy velocity of  m/s) of the reservoir was successfully modeled, in agreement with what was inferred by the extensive deep explorations campaigns performed in the area.

Simulation of the exploitation process covered a total time interval of 10,000 years with fluid extraction and injection limited to the initial 50 years. Simulations showed that only small changes were induced by the exploitation of the geothermal system (producing well temperature increase between 2.5 and 9.5°C after 50 years) and no thermal breakthrough occurs. Full recovery occurs in about one thousand years due to the highly convective behavior of the reservoir. The good agreement between measured and simulated results for the natural state allowed a confident prediction of the reservoir response to future exploitation.

These concluding remarks were also sustained by the qualitatively similar calculated patterns resulting from the FEFLOW performed simulation. Even if a time discrepancy in the identification of the present-day natural state occurs between FEFLOW and OGS, the convective system behavior, the fitting between simulated and real thermal data, and the reservoir response to the tested exploitation scenario are fully comparable.

Such models support the understanding of reservoir behavior and are critical to optimal reservoir management and sustainable utilization. Their reliability could be improved by integrating data from new superficial and deep explorations, and, at the same time, they can support the planning of new investigations, well drilling, and the design of exploitation steps aimed at the usage of geothermal energy in the Caste Giorgio-Torre Alfina area.


This work represents the continuation of a project that was introduced in its early stage with the oral presentation at 9th EGU General Assembly (Geophysical Research Abstracts Vol. 19, EGU2017-16968).

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This research originated from a collaboration between University of Milano-Bicocca and Ricerca di Sistema Energetico (RSE). The authors acknowledge RSE for providing the data, in particular, Fabio Moia who encouraged this cooperation. They also greatly acknowledge Roberto Guandalini and Agate Giordano, from the RSE research group, for discussing with them the modelling approach, and for the helpful suggestions in the construction and the calibration of the model.