Research Article  Open Access
Giorgio Volpi, Fabien Magri, Francesca Colucci, Thomas Fisher, Mattia De Caro, Giovanni B. Crosta, "Modeling Highly Buoyant Flows in the Castel Giorgio: Torre Alfina Deep Geothermal Reservoir", Geofluids, vol. 2018, Article ID 3818629, 19 pages, 2018. https://doi.org/10.1155/2018/3818629
Modeling Highly Buoyant Flows in the Castel Giorgio: Torre Alfina Deep Geothermal Reservoir
Abstract
The Castel GiorgioTorre Alfina (CGTA, 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 CGTA to simulate the undisturbed natural geothermal field and investigate the impacts of the exploitation process. The open source finiteelement 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 50year 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 CGTA 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, windenergy, and several forms of solar energy, such as bioenergy, biofuel, photovoltaic, and solarthermal 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 [4–6], northern California [7], the Pannonic Basin [8], and the Rhine Graben [9–11]. 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, 12–18].
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 heatflux 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/m^{2} (Larderello system) and 600 mW/m^{2} (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 GiorgioTorre Alfina field (CGTA, northern Latium, Figure 1) is an example of an early explored and so far not exploited mediumenthalpy geothermal system [23–25]. Detailed hydrogeothermal data, available for the selected area since early 70s, show that the CGTA is a potential geothermal reservoir with medium thermal characteristics (120°C–210°C) whose fluids (pressurized water and gas, mainly CO_{2}) are hosted in a fractured carbonate formation [24, 26–30]. 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 CGTA area suitable for future exploitation through a new generation 5 MWe geothermal pilot power plant. Following the guidelines of the abovementioned 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, 32–43]. 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, mediumenthalpy CGTA reservoir to reproduce the highly convective undisturbed presentday 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 50year operational life cycle adopting a well doublet system at a 1,050 t/h flow rate [48]. The finiteelement 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 finiteelement 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 CGTA reservoir. Then, the numerical approach based on the OpenGeoSys software is given. Results are obtained both at shortterm (i.e., operational) and longterm (i.e., full reservoir recovery) time scale. Besides providing valuable guidelines for future exploitation of the CGTA deep geothermal reservoir, this study highlights the importance of field data constraints for the interpretation of numerical results of fluid processes in reservoirscale systems.
2. Reservoir Characterization
2.1. Regional Geological Setting
The occurrence of medium and highenthalpy 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 ([51–54] and references therein).
The presentday 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 AlpineApennine orogenesis [18, 55]. The compressive phase resulted in the formation of foldandthrustbelts and associated piggyback basins with NNESSW oriented trend [56–58]. Then, the subsequent extensional phase, due to the Tyrrhenian backarc 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/subLigurian complex: JurassicEocene clayeymarly 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: TriassicLower Miocene arenaceous and clayeymarly formations, calcareoussiliceous 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 GiorgioTorre Alfina Geothermal Field
The CGTA 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 CGTA 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 NS 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 CGTA field (19711972) was aimed to reach and cross the argillaceous and shaly terrains of the Ligurian and SubLigurian 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 CO_{2}. This 100meter thick cap, recognized only in the central part of the field, was extensively exploited until few years ago for CO_{2} storage by the well Alfina013 (Figure 1).
The target of the more recent campaign (19871988) 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 knickpoint 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, largescale convective flow confined in the area of the buried structural high.
In summary, the CGTA 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 km^{2} (Figure 4), is located north of the Vulsini calderas and it is bordered by the MesoCenozoic 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 EW direction and of 2.5 km in the NS 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 × 20meter 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 flyschtype 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 strikeslip and subordinate normal fault systems (with associated fracture network) cutting and dislocating the internal architecture of the reservoir [60]. No anomalous soil CO_{2} 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/subLigurian 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 km^{2}) and those falling outside the producing area. The model includes also a NESW 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 reservoirscale threedimensional thermohydraulic (TH) model was built to investigate the different processes involved in the CGTA geothermal reservoir.
The open source finiteelement simulator OpenGeoSys (OGS) [49] was used to solve the differential equations governing densitydriven 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 presentday, highly convective, unexploited, natural state of the CGTA 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 longterm (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 finiteelement 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 heatflux of 0.256 W/m^{2} (i.e., Neumann type) was applied at the nodes on bottom boundary below the reservoir area (Figure 5).
A nomass 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 presentday reservoir condition, without any fluid extraction/injection scenarios. To let the system reach the presentday 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 50year 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 timedependent. At the end of the 50year 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 heatflux 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 [23–25, 30, 31, 57, 61–63, 65, 70]. Default values for thermal conductivity of water (0.65 W/m/K) and heat capacity of water (4.2 MJ/m^{3}/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.
(a) 0 years
(b) 20,000 years
(c) 125,000 years
4. Natural State Simulation
4.1. Model Definition
The natural state simulation aimed to define the presentday, 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 m^{2}, 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 }m^{2} (Table 1). A very low permeability value of 10^{–18 }m^{2} 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 upflows, with a fluid velocity in the range of 2–4 10^{8 }m/s, and associated lateral downflows. 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 bestfitting 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,000year 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 welldeveloped 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 temperaturedepth 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 bestfitting 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 CGTA 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 presentday distribution of both temperature and pressure field, were imported from the calculated bestfitting 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 longterm 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).
(a)
(b)
(c)
(d)
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.
(a)
(b)
(c)
(d)
(e)
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 CGTA reservoir to the production process, the influence area of the “coldwater” 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 CGTA reservoir and then simulate the exploitation of this undeveloped geothermal field. A general procedure for model calibration was applied [25, 35, 73–75], 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 bestfitting simulation time (125,000 years, Figure 11), the OGS model was sliced with a vertical plane along the EW axis and passing through the Alfina015 well (AA′ 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 crosssection plane at different relative positions with respect to the convective plumes (i.e., interplume or inbetween 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 convectiondominated geothermal systems, for which comparable analyses were performed [18, 22, 43, 47, 76, 77].
Multiple horizontal temperature profiles along AA′ (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 presentday 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 longterm 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., bestfitting 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 bestfitting 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 presentday situation modeled with the two software programs (i.e., bestfitting 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 bestfitting 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 presentday 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 presentday 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)).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Finally, the areal extent of the 80°C “coldwater” 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 CGTA mediumenthalpy geothermal field, to investigate the feasibility of geothermal exploitation, and to test capabilities of different codes at modeling highly buoyant flows. A fitforpurpose 3D numerical model of the CGTA 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 presentday, 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 presentday 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 GiorgioTorre Alfina area.
Disclosure
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, EGU201716968).
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research originated from a collaboration between University of MilanoBicocca 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.
References
 L. Rybach, “Ingrid Stober und Kurt Bucher: Geothermie,” uwf UmweltWirtschaftsForum, vol. 20, no. 24, pp. 197–197, 2012. View at: Publisher Site  Google Scholar
 M. H. Dickson and M. Fanelli, “Small Geothermal Resources: A Review,” Energy Sources, vol. 16, no. 3, pp. 349–376, 2007. View at: Publisher Site  Google Scholar
 International Energy Agency, “Energy statistics of NonOECD countries,” Tech. rep, 2009. View at: Google Scholar
 G. P. Eaton, R. L. Christiansen, H. M. Iyer et al., “Magma beneath Yellowstone National Park,” Science, vol. 188, no. 4190, pp. 787–796, 1975. View at: Publisher Site  Google Scholar
 P. Morgan, D. D. Blackwell, R. E. Spafford, and R. B. Smith, “Heat flow measurements in Yellowstone Lake and the thermal structure of the Yellowstone Caldera,” Journal of Geophysical Research: Atmospheres, vol. 82, no. 26, pp. 3719–3732, 1977. View at: Publisher Site  Google Scholar
 I. Lucchitta, “Role of heat and detachment in continental extension as viewed from the eastern basin and range province in Arizona,” Tectonophysics, vol. 174, no. 12, pp. 77–114, 1990. View at: Publisher Site  Google Scholar
 L. W. Younker, P. W. Kasameyer, and J. D. Tewhey, “Geological, geophysical, and thermal characteristics of the Salton Sea Geothermal Field, California,” Journal of Volcanology and Geothermal Research, vol. 12, no. 34, pp. 221–258, 1982. View at: Publisher Site  Google Scholar
 D. Ravnik, D. Rajver, M. Poljak, and M. Živčić, “Overview of the geothermal field of Slovenia in the area between the Alps, the Dinarides and the Pannonian basin,” Tectonophysics, vol. 250, no. 13, pp. 135–149, 1995. View at: Publisher Site  Google Scholar
 D. Werner and H. Kahle, “A geophysical study of the Rhinegraben — I. Kinematics and geothermics,” The Geophysical Journal of the Royal Astronomical Society, vol. 62, no. 3, pp. 617–629, 1980. View at: Publisher Site  Google Scholar
 J. P. Brun, M.A. Gutscher, and dekorpecors teams, “Deep crustal structure of the Rhine Graben from dekorpecors seismic reflection data: A summary,” Tectonophysics, vol. 208, no. 13, pp. 139–147, 1992. View at: Publisher Site  Google Scholar
 S. Bellani, A. Brogi, A. Lazzarotto, D. Liotta, and G. Ranalli, “Heat flow, deep temperatures and extensional structures in the Larderello Geothermal Field (Italy): Constraints on geothermal fluid flow,” Journal of Volcanology and Geothermal Research, vol. 132, no. 1, pp. 15–29, 2004. View at: Publisher Site  Google Scholar
 A. Barelli, R. Bertani, G. Cappetti, and A. Ceccarelli, “An update on Travale Radicondoli geothermal field,” in Proceedings of the Proceedings World Geothermal Congress, pp. 1581–1586, Florence , Italy, 1995. View at: Google Scholar
 A. Barelli, G. Cappetti, and G. Stefani, “Optimum exploitation strategy at LarderelloValle Secolo,” in Proceedings of the Proceedings World Geothermal Congress, pp. 1779–1783, Florence , Italy, 1995. View at: Google Scholar
 A. Barelli, G. Cappetti, and G. Stefani, “Results of deep drilling in the LarderelloTravale/Radicondoli geothermal area,” in Proceedings of the Proceedings World Geothermal Congress, pp. 1275–1278, Florence , Italy, 1995. View at: Google Scholar
 A. Barelli, G. Bertini, G. Buonasorte, G. Cappetti, and A. Fiordelisi, “Recent deep exploration results at the margins of the LarderelloTravale geothermal system,” in Proceedings of the Proceedings World Geothermal Congress, pp. 965–970, KyushuTohoku, Japan, 2000. View at: Google Scholar
 F. Batini, A. Brogi, A. Lazzarotto, D. Liotta, and E. Pandeli, “Geological features of LarderelloTravale and Mt. Amiata geothermal areas (southern Tuscany, Italy),” Episodes, vol. 26, no. 3, pp. 239–244, 2003. View at: Google Scholar
 G. N. Tiwari and M. K. Ghosal, “Renewable energy resources: Basic principles and applications,” in Proceedings of the International Journal of Industrial Engineering Computations, vol. 3, pp. 649–662, 2005. View at: Google Scholar
 P. Romagnoli, A. Arias, A. Barelli, M. Cei, and M. Casini, “An updated numerical model of the LarderelloTravale geothermal system, Italy,” Geothermics, vol. 39, no. 4, pp. 292–313, 2010. View at: Publisher Site  Google Scholar
 L. Carmignani, F. A. Decandia, P. L. Fantozzi, A. Lazzarotto, D. Liotta, and M. Meccheri, “Tertiary extensional tectonics in Tuscany (Northern Apennines, Italy),” Tectonophysics, vol. 238, no. 14, pp. 295–315, 1994. View at: Publisher Site  Google Scholar
 C. Brunet, P. Monié, L. Jolivet, and J.P. Cadet, “Migration of compression and extension in the Tyrrhenian Sea, insights from 40Ar/39Ar ages on micas along a transect from Corsica to Tuscany,” Tectonophysics, vol. 321, no. 1, pp. 127–155, 2000. View at: Publisher Site  Google Scholar
 B. Della Vedova, S. Bellani, G. Pellis, and P. Squarci, “Deep temperatures and surface heat flow distribution,” in Anatomy of an orogen: the Apennines and adjacent Mediterranean basins, pp. 65–76, 2001. View at: Google Scholar
 B. Della Vedova, C. Vecellio, S. Bellani, and U. Tinivella, “Thermal modelling of the Larderello geothermal field (Tuscany, Italy),” International Journal of Earth Sciences, vol. 97, no. 2, pp. 317–332, 2008. View at: Publisher Site  Google Scholar
 R. Cataldi and M. Rendina, “Recent discovery of a new geothermal field in Italy: Alfina,” Geothermics, vol. 2, no. 34, pp. 106–116, 1973. View at: Publisher Site  Google Scholar
 G. Buonasorte, R. Cataldi, A. Ceccarelli et al., “Ricerca ed esplorazione nell'area geotermica di Torre Alfina (LazioUmbria),” in Bollettino della Società Geologica Italiana, vol. 107, pp. 265–337, 1988. View at: Google Scholar
 Colucci and Guandalini, “Modelli geologici e simulazione numerica di sistemi geotermici,” Rapporto di Ricerca di Sistema, 2014, http://doc.rseweb.it/doc/docsfoglia/15000985316054/15000985316054.html. View at: Google Scholar
 F. Barberi, G. Buonasorte, R. Cioni et al., “PlioPleistocene geological evolution of the geothermal area of Tuscany and Latium,” Mem. Descr. Carta Geol. Ital, vol. 49, pp. 77–133, 1994. View at: Google Scholar
 C. Chiarabba, A. Amato, and A. Fiordelisi, “Upper crustal tomographic images of the AmiataVulsini geothermal region, central Italy,” Journal of Geophysical Research: Atmospheres, vol. 100, no. 3, pp. 4053–4066, 1995. View at: Publisher Site  Google Scholar
 G. Chiodini, F. Frondini, and F. Ponziani, “Deep structures and carbon dioxide degassing in Central Italy,” Geothermics, vol. 24, no. 1, pp. 81–94, 1995. View at: Publisher Site  Google Scholar
 M. Doveri, M. Lelli, L. Marini, and B. Raco, “Revision, calibration, and application of the volume method to evaluate the geothermal potential of some recent volcanic areas of Latium, Italy,” Geothermics, vol. 39, no. 3, pp. 260–269, 2010. View at: Publisher Site  Google Scholar
 M. L. Carapezza, M. Ranaldi, A. Gattuso, N. M. Pagliuca, and L. Tarchini, “The sealing capacity of the cap rock above the Torre Alfina geothermal reservoir (Central Italy) revealed by soil CO2 flux investigations,” Journal of Volcanology and Geothermal Research, vol. 291, pp. 25–34, 2015. View at: Publisher Site  Google Scholar
 G. Buonasorte, E. Pandeli, and A. Fiordelisi, “The Alfina 15 Well; deep geological data from northern Latium (Torre Alfina geothermal area),” Bollettino della Societa Geologica Italiana, pp. 823–831, 1991. View at: Google Scholar
 E. U. Antunez, S. K. Sanyal, A. J. Menzies et al., “Forecasting well and reservoir behavior using numerical simulation, Uenotai geothermal field, Akita prefecture, Japan,” in Proceedings of the 1990 International Symposium on Geothermal Energy, pp. 1255–1262, August 1990. View at: Google Scholar
 E. U. Antúnez, G. S. Bodvarsson, and M. A. Walters, “Numerical simulation study of the Northwest Geysers geothermal field, a case study of the Coldwater Creek steamfield,” Geothermics, vol. 23, no. 2, pp. 127–141, 1994. View at: Publisher Site  Google Scholar
 M. O'Sullivan, B. Barnett, and Y. Razali, “Numerical simulation of the Kamojang Geothermal Field, Indonesia,” in Proceedings of the 1990 International Symposium on Geothermal Energy, pp. 1317–1324, August 1990. View at: Google Scholar
 M. J. O'Sullivan, K. Pruess, and M. J. Lippmann, “State of the art geothermal reservoir simulation,” Geothermics, vol. 30, no. 4, pp. 395–429, 2001. View at: Publisher Site  Google Scholar
 M. Hanano, “Reservoir engineering studies of the matsukawa geothermal field, Japan,” in Proceedings of the 1992 Annual Meeting of the Geothermal Resources Council, pp. 643–650, October 1992. View at: Google Scholar
 M. Hanano, “Simulation study of the Matsukawa geothermal reservoir: Natural state and its response to exploitation,” Journal of Energy Resources TechnologyTransactions of the ASME, vol. 114, no. 4, pp. 309–314, 1992. View at: Publisher Site  Google Scholar
 G. Axelsson and G. Bjornsson, “Detailed threedimensional modeling of the Botn hydrothermal system in NIceland,” in Proceedings of the 18th Workshop on Geothermal Reservoir Engineering, pp. 159–166, Stanford University, Stanford, California, 1993. View at: Google Scholar
 M. Pham and A. J. Menzies, “Results from a fieldwide numerical model of the geysers geothermal field, California,” in Proceedings of the 1993 Annual Meeting on Utilities and Geothermal: An Emerging Partnership, pp. 259–265, October 1993. View at: Google Scholar
 M. Pham, A. J. Menzies, S. K. Sanyal et al., “Numerical modeling of the hightemperature geothermal system of Amatitlan, Guatemala,” in Proceedings of the 1996 Annual Meeting of the Geothermal Resources Council, pp. 833–838, October 1996. View at: Google Scholar
 R. Bertani and G. Cappetti, “Numerical simulation of the Monteverdi zone (western border of the Larderello geothermal field,” in Proceedings of the Proceedings World Geothermal Congress 95, pp. 1735–1740, Florence, 1995. View at: Google Scholar
 M. A. Antics, “Computer simulation of the Oradea geothermal reservoir,” in Proceedings of the 22nd Workshop on Geothermal Reservoir Engineering, pp. 491–495, Stanford University, Stanford, Calif, USA, 1997. View at: Google Scholar
 P. Fulignati, P. Marianelli, A. Sbrana, and V. Ciani, “3D geothermal modelling of the mount amiata hydrothermal system in Italy,” Energies, vol. 7, no. 11, pp. 7434–7453, 2014. View at: Publisher Site  Google Scholar
 T. Li, S. Shiozawa, and M. W. McClure, “Thermal breakthrough calculations to optimize design of a multiplestage Enhanced Geothermal System,” Geothermics, vol. 64, pp. 455–465, 2016. View at: Publisher Site  Google Scholar
 M. Caputo, “Diffusion of fluids in porous media with memory,” Geothermics, vol. 28, no. 1, pp. 113–130, 1999. View at: Publisher Site  Google Scholar
 K. Nicholson, Geothermal fluids: chemistry and exploration techniques, Springer Science & Business Media, 2012.
 M. G. Blöcher, G. Zimmermann, I. Moeck, W. Brandt, A. Hassanzadegan, and F. Magri, “3D numerical modeling of hydrothermal processes during the lifetime of a deep geothermal reservoir,” Geofluids, vol. 10, no. 3, pp. 406–421, 2010. View at: Publisher Site  Google Scholar
 TW & LKW Geotermia Italia S.p.A., Impianto pilota geotermico Castel Giorgio (TR). Progetto definitivo e programma lavori, 2013, http://www.va.minambiente.it/itIT/Oggetti/Documentazione/1373/1855?Testo=%20Progetto%20e%20Programma%20Lavori&pagina=1#formcercaDocumentazione.
 O. Kolditz, S. Bauer, L. Bilke et al., “OpenGeoSys: An opensource initiative for numerical simulation of thermohydromechanical/chemical (THM/C) processes in porous media,” Environmental Earth Sciences, vol. 67, no. 2, pp. 589–599, 2012. View at: Publisher Site  Google Scholar
 H.J. G. Diersch, Feflow: Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media, SpringerVerlag Berlin Heidelberg, Berlin, Germany, 2014. View at: Publisher Site
 F. A. Decandia, A. Lazzarotto, D. Liotta, L. Cernobori, and R. Nicolich, “traverse: insights on postcollisional evolution of Northern Apennines,” The CROP, vol. 03, pp. 427–439, 1998. View at: Google Scholar
 D. Liotta, L. Cernobori, and R. Nicolicl, “Restricted rifting and its coexistence with compressional structures: Results from the CROP 3 traverse (northern Apennines, Italy),” Terra Nova, vol. 10, no. 1, pp. 16–20, 1998. View at: Publisher Site  Google Scholar
 A. Brogi, A. Lazzarotto, D. Liotta, and G. Ranalli, “Extensional shear zones as imaged by reflection seismic lines: The Lardello geothermal field (central Italy),” Tectonophysics, vol. 363, no. 12, pp. 127–139, 2003. View at: Publisher Site  Google Scholar
 A. Brogi, A. Lazzarotto, D. Liotta et al., “Crustal structures in the geothermal areas of southern Tuscany (Italy): Insights from the CROP 18 deep seismic reflection lines,” Journal of Volcanology and Geothermal Research, vol. 148, no. 12, pp. 60–80, 2005. View at: Publisher Site  Google Scholar
 E. Carminati and C. Doglioni, “Mediterranean tectonics,” Encyclopedia of Geology, vol. 2, pp. 135–146, 2005. View at: Google Scholar
 G. Buonasorte, A. Fiordelisi, E. Pandeli, U. Rossi, and F. Sollevanti, “Stratigraphic correlations and structural setting of the pre neoautochthonous sedimentary sequences of northern Latium,” Periodico di Mineralogia, vol. 56, no. 23, pp. 111–122, 1987. View at: Google Scholar
 G. Buonasorte, G. M. Cameli, A. Fiordelisi, M. Parotto, and I. Perticone, “Results of geothermal exploration in Central Italy (LatiumCampania),” in Proceedings of the In Proceedings of the World Geothermal Congress, pp. 18–31, Florence, Italy, 1995. View at: Google Scholar
 L. Petracchini, D. Scrocca, S. Spagnesi, and F. Minelli, “3D geological modeling to support the assessment of conventional and unconventional geothermal resources in the latium region (Central Italy),” in Proceedings of the In World Geothermal Congress, pp. 19–25, 2015. View at: Google Scholar
 A. Minissale, “The Larderello geothermal field: a review,” EarthScience Reviews, vol. 31, no. 2, pp. 133–151, 1991. View at: Publisher Site  Google Scholar
 F. Rossetti, C. Faccenna, L. Jolivet, R. Funiciello, F. Tecce, and C. Brunet, “Syn versus postorogenic extension: The case study of Giglio Island (Northern Tyrrhenian Sea, Italy),” Tectonophysics, vol. 304, no. 12, pp. 71–93, 1999. View at: Publisher Site  Google Scholar
 L. D. Nardi, G. Pieretti, and M. Rendina, “Stratigrafia dei terreni perforati dai sondaggi ENEL nell9area geotermica di Torre Alfina,” Bollettino della Societa Geologica Italiana, vol. 96, no. 3, pp. 403–422, 1977. View at: Google Scholar
 D. Cosentino, P. Cipollari, P. Marsili, and D. Scrocca, “Geology of the central Apennines: A regional review,” Journal of the Virtual Explorer, vol. 36, 2010. View at: Publisher Site  Google Scholar
 G. Vignaroli, A. Pinton, A. A. De Benedetti et al., “Structural compartmentalisation of a geothermal system, the Torre Alfina field (central Italy),” Tectonophysics, vol. 608, pp. 482–498, 2013. View at: Publisher Site  Google Scholar
 D. Piscopo, M. Gattiglio, E. Sacchi, and E. Destefanis, “Tectonicallyrelated fluid circulation in the san casciano dei bagnisarteano area (m. cetona ridgesouthern tuscany): a coupled structural and geochemical investigation,” Bollettino della Societa Geologica Italiana, vol. 128, no. 2, pp. 575–586, 2009. View at: Publisher Site  Google Scholar
 A. Costantini, C. Ghezzo, and A. Lazzarotto, “Carta geologica dell'area geotermica di Torre Alfina (prov. di SienaViterboTerni). Ente Nazionale per l'Energia Elettrica (ENEL, Unitα Nazionale Geotermica Pisa. Cartografia S.EL.C.A,” Carta geologica dell'area geotermica di Torre Alfina (prov. di SienaViterboTerni). Ente Nazionale per l'Energia Elettrica (ENEL, Unitα Nazionale Geotermica Pisa. Cartografia S.EL.C.A, 1984. View at: Google Scholar
 L. Carmignani and A. Lazzarotto, “Carta geologica della Toscana / geological Map of Tuscany (Italy,” in Regione Toscana, Direzione delle Politiche Territoriali e AmbientaliServizio geologico, vol. 1, Carta geologica della Toscana / geological Map of Tuscany (Italy) 1, 250000. Regione Toscana, 2004. View at: Google Scholar
 F. Magri, N. Inbar, C. Siebert, E. Rosenthal, J. Guttman, and P. Möller, “Transient simulations of largescale hydrogeological processes causing temperature and salinity anomalies in the Tiberias Basin,” Journal of Hydrology, vol. 520, pp. 342–355, 2015. View at: Publisher Site  Google Scholar
 EMSI, Groundwater Modeling System 6.0. Environmental Modeling Systems, Inc, http://www.aquaveo.com/software/gmsgroundwatermodelingsystemintroduction, 2006.
 L. Marini, F. Franceschini, M. Ghigliotti, M. Guidi, and A. Merla, “Valutazione del potenziale geotermico nazionale. ENEAGeotermica Italiana Report for the Ministero dellIndustria,” del Commercio e dellArtigianato, 1993. View at: Google Scholar
 G. Giordano, A. A. De Benedetti, A. Diana et al., “The Colli Albani mafic caldera (Roma, Italy): Stratigraphy, structure and petrology,” Journal of Volcanology and Geothermal Research, vol. 155, no. 12, pp. 49–80, 2006. View at: Publisher Site  Google Scholar
 E. A. Porras, T. Tanaka, H. Fujii, and R. Itoi, “Numerical modeling of the Momotombo geothermal system, Nicaragua,” Geothermics, vol. 36, no. 4, pp. 304–329, 2007. View at: Publisher Site  Google Scholar
 E. M. Llanos, S. J. Zarrouk, and R. A. Hogarth, “Numerical model of the Habanero geothermal reservoir, Australia,” Geothermics, vol. 53, pp. 308–319, 2015. View at: Publisher Site  Google Scholar
 M. J. O'sullivan, “Geothermal reservoir simulation,” International Journal of Energy Research, vol. 9, no. 3, pp. 319–332, 1985. View at: Publisher Site  Google Scholar
 G. S. Bodvarsson, K. Pruess, V. Stefansson et al., Natural state model of the Nesjavellir geothermal field, Iceland (No. SGPTR9317). Lawrence Berkeley National Laboratory (LBNL), Earth Sciences Division, Berkeley, CA, USA, National Energy Authority of Iceland, Reykjavik, Iceland; Reykjavik Municipal District Heating Service, Reykjavik, Iceland, 1986.
 K. Pruess, J. S. Y. Wang, and Y. W. Tsang, “On thermohydrologic conditions near high‐level nuclear wastes emplaced in partially saturated fractured tuff: 1. Simulation studies with explicit consideration of fracture effects,” Water Resources Research, vol. 26, no. 6, pp. 1235–1248, 1990. View at: Publisher Site  Google Scholar
 B. M. Feather and R. C. M. Malate, “Numerical modeling of the Mita geothermal field, Cerro Blanco, Guatemala,” in Proceedings of the In Proceedings of the ThirtyEighth Workshop on Geothermal Reservoir Engineering, pp. 11–13, Stanford University, Stanford, CA, February 2013. View at: Google Scholar
 A. Ebigbo, J. Niederau, G. Marquart et al., “Influence of depth, temperature, and structure of a crustal heat source on the geothermal reservoirs of Tuscany: numerical modelling and sensitivity study,” Geothermal Energy, vol. 4, no. 1, article no. 5, 2016. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Giorgio Volpi 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.