The Integration of 3 D Modeling and Simulation to Determine the Energy Potential of Low-Temperature Geothermal Systems in the Pisa ( Italy ) Sedimentary Plain

Shallow, low-temperature geothermal resources can significantly reduce the environmental impact of heating and cooling. Based on a replicable standard workflow for three-dimensional (3D) geothermal modeling, an approach to the assessment of geothermal energy potential is proposed and applied to the young sedimentary basin of Pisa (north Tuscany, Italy), starting from the development of a geothermal geodatabase, with collated geological, stratigraphic, hydrogeological, geophysical and thermal data. The contents of the spatial database are integrated and processed using software for geological and geothermal modeling. The models are calibrated using borehole data. Model outputs are visualized as three-dimensional reconstructions of the subsoil units, their volumes and depths, the hydrogeological framework, and the distribution of subsoil temperatures and geothermal properties. The resulting deep knowledge of subsoil geology would facilitate the deployment of geothermal heat pump technology, site selection for well doublets (for open-loop systems), or vertical heat exchangers (for closed-loop systems). The reconstructed geological–hydrogeological models and the geothermal numerical simulations performed help to define the limits of sustainable utilization of an area’s geothermal potential.


Introduction
Sustainability, including in urban areas, has become a key concern when planning for the future [1].Increasing energy demand, the changing climate, and energy resources limitations have prompted several countries to introduce strategies to enhance energy-efficiency and production sustainability.The construction and ongoing use of the built environment are substantial contributors to growing energy demand, using from 20 to 40% of the energy consumed in developed countries and accounting for more than 30% of total greenhouse emissions [1,2].In Europe, the Directive on Energy Performance of Buildings [3] states that all new constructions should be "zero-energy buildings" by the end of 2020.This implies buildings that consume minimum levels of renewably sourced energy.By replacing fossil fuel sources and substantially reducing the CO 2 emissions, renewable energy sources help interpolation algorithm consists of three steps: (1) Refine-increase grid resolution; (2) Snap-regrid the data; and (3) Smooth-minimize curvature of the grid [18].
Starting with the upper surfaces, a 3D grid was used as the basis of a 3D geological model of the investigated area.The initial conceptual model was revised using improved 3D elaboration and visualization: For example, lithostratigraphic and well log data can suggest adjustments to the geothermal geodatabase, while the 3D reconstruction of aquifers allows the conceptual model to be further modified and improved.The refined geological conceptual model and the resultant 3D Petrel subsurface model became the framework for subsequent geothermal simulations.The top surfaces of the subsurface units, from the 3D modeling, were imported into TOUGH2 [19] through the PetraSim (version 5; Rockware Inc, Golden, CO, USA) graphical interface.TOUGH2 is a non-isothermal multiphase fluid flow simulator for modeling heat and mass transfer in geothermal systems.In this study, the EOS1 module was used.The computational approach used in TOUGH2 is based on the integral finite difference (IFD) method [20] using fully implicit time discretization.The residual form of the system of equations, considering the underground lithological setting, the heat flow dynamics, and the heat exchange and fluid mass transport laws, is solved using the Newton-Raphson iteration method [21].Geothermal parameters, including density, porosity, permeability, heat conductivity, and specific heat (both literature-derived values and in situ measurements), were assigned as input data to subsoil units.The boundary conditions were set as for a closed system.A grid size was chosen for the calculation based on the required spatial resolution and the acceptable computational load.The simulation then modeled the steady state of the subsoil volume.
Temperature data, as measured in water, geothermal, and deep oil and gas exploratory wells, were used to check the thermal settings of the model.Several iterations were performed to obtain the best fit between measured and simulated temperature data.The distribution of output temperatures, obtained through the TOUGH2 simulation, was transferred back to the Petrel software, allowing visualization of the 3D distribution of subsoil temperatures in the aquifer or geological units (along the X-, Y-, or Z-axis, in 2D sections or as well logs).Further, the simulated temperatures may be integrated with temperature values measured in existing wells, giving a robust picture of the 3D temperature distribution in the subsoil.These temperatures were interpolated to a continuous raster surface, using an inverse distance weighted (IDW) interpolation [22].IDW is a deterministic interpolation method that give greater weight to points closer to the unknown point that we want to estimate, than to distant ones.Cell values are averaged over the nearest data points, and the influence of distant data decreases exponentially, with the power of 2. This process allows the production of temperature maps at different depths.
Moreover, the GIS raster calculator tool was used to generate thematic maps.It allows execution of map algebra expressions involving two or more rasters, resulting in a new raster.In this study, starting from the top surface of the lithological units (input raster) and using the "volumetric method equation" [23,24], the volumes of the lithological units were determined.The outputs allow easy estimation and mapping of geothermal potential and other useful geothermal parameters, including temperatures and aquifer volumes.Next, a conceptual model of the subsoil was defined for the area, based on the encoded lithostratigraphic, hydrogeological and geothermal knowledge, together with the scientific and technical literature.
Subsurface modeling is normally based on interpretation of two-dimensional (2D) cross-sections, maps, well logs, and geophysical surveys.The incorporation of 3D modeling allows interpolation of the subsurface units and their features, significantly improving understanding and visualization of the subsoil geological framework.In our approach, the integration of point data, 2D cross-sections and 3D modeling was performed using the Schlumberger Petrel E and P Software Platform (v.2016.1,Schlumberger-Software Integrated Solutions, Parma, Italy).This software is a complete application for 3D geological-geophysical modeling.
All required data were imported from the designed GIS geodatabase.Point data were converted into areal information using the Petrel "convergent interpolation algorithm", an iterative algorithm that delineates the top surfaces of the main lithological units.Each iteration of the convergent interpolation algorithm consists of three steps: (1) Refine-increase grid resolution; (2) Snap-regrid the data; and (3) Smooth-minimize curvature of the grid [18].
Starting with the upper surfaces, a 3D grid was used as the basis of a 3D geological model of the investigated area.The initial conceptual model was revised using improved 3D elaboration and visualization: For example, lithostratigraphic and well log data can suggest adjustments to the geothermal geodatabase, while the 3D reconstruction of aquifers allows the conceptual model to be further modified and improved.The refined geological conceptual model and the resultant 3D Petrel subsurface model became the framework for subsequent geothermal simulations.The top surfaces of the subsurface units, from the 3D modeling, were imported into TOUGH2 [19] through the PetraSim (version 5; Rockware Inc, Golden, CO, USA) graphical interface.TOUGH2 is a non-isothermal multiphase fluid flow simulator for modeling heat and mass transfer in geothermal systems.In this study, the EOS1 module was used.The computational approach used in TOUGH2 is based on the integral finite difference (IFD) method [20] using fully implicit time discretization.The residual form of the system of equations, considering the underground lithological setting, the heat flow dynamics, and the heat exchange and fluid mass transport laws, is solved using the Newton-Raphson iteration method [21].Geothermal parameters, including density, porosity, permeability, heat conductivity, and specific heat (both literature-derived values and in situ measurements), were assigned as input data to subsoil units.The boundary conditions were set as for a closed system.A grid size was chosen for the calculation based on the required spatial resolution and the acceptable computational load.The simulation then modeled the steady state of the subsoil volume.
Temperature data, as measured in water, geothermal, and deep oil and gas exploratory wells, were used to check the thermal settings of the model.Several iterations were performed to obtain the best fit between measured and simulated temperature data.The distribution of output temperatures, obtained through the TOUGH2 simulation, was transferred back to the Petrel software, allowing visualization of the 3D distribution of subsoil temperatures in the aquifer or geological units (along the X-, Y-, or Z-axis, in 2D sections or as well logs).Further, the simulated temperatures may be integrated with temperature values measured in existing wells, giving a robust picture of the 3D temperature distribution in the subsoil.These temperatures were interpolated to a continuous raster surface, using an inverse distance weighted (IDW) interpolation [22].IDW is a deterministic interpolation method that give greater weight to points closer to the unknown point that we want to estimate, than to distant ones.Cell values are averaged over the nearest data points, and the influence of distant data decreases exponentially, with the power of 2. This process allows the production of temperature maps at different depths.
Moreover, the GIS raster calculator tool was used to generate thematic maps.It allows execution of map algebra expressions involving two or more rasters, resulting in a new raster.In this study, starting from the top surface of the lithological units (input raster) and using the "volumetric method equation" [23,24], the volumes of the lithological units were determined.The outputs allow easy estimation and mapping of geothermal potential and other useful geothermal parameters, including temperatures and aquifer volumes.

Case History
The Pisa plain consists primarily of neoautochthonous deposits that fill a wide and deep half-graben striking NW-SE [25].The Pleistocene evolution of the sedimentary basin has resulted in a complex stratigraphic pattern in which marine and continental sand-gravel layers alternate and interdigitate with clay layers [26] covering a bedrock formed by carbonate units pertaining to the Tuscan Nappe [27].

Database Setting, Data Screening, Collection, and Processing
The available data, including maps and databases of existing wells with subsoil temperatures and water composition, for the Pisa sedimentary plain were collected from diverse sources.The data were often not in digital form and varied in nature and quality.To define a trustworthy conceptual geothermal model of the area, we undertook a new hydrogeochemical survey that provided accurate contemporary information about aquifer quality and temperature distributions in the subsoil.
Direct temperature measurement in existing water wells is particularly difficult because of physical obstacles (such as well heads or submersible pumps that make it impossible to insert thermometers); however, these data are fundamental to the subsequent steps of the workflow (the thermal simulation being calibrated on measured temperatures).To resolve the problem, we selected water wells having a constrained depth of water captation (filters or slotted liners) in the multilayer aquifer.The temperature of the aquifer was recorded at the well head after the observed water temperature stabilized.This is valid, as the deviation between temperatures measured at the well head, in dynamic conditions, and direct measurements in the aquifer has been found to be less than 1%.
The data storage process requires that the data are properly structured and encoded.A geothermal database was designed, organized, and implemented, and four information layers were mutually connected using a key identity field.Geographic coordinates were converted into a common reference projection, WGS_1984_UTM_Zone_32N (EPSG 32632), in the GIS [41].The GIS also contained raster elements including geological maps, thematic maps and a digital terrain model (DTM).

Geological Conceptual Model
The lithostratigraphy of the shallower sedimentary sequence, to a depth of about 200 m, was reconstructed using the borehole database (Figure 2a).This corresponded to the maximum depth reached by water wells in the area; higher depths are reached only by a few other wells (geothermal and hydrocarbon exploration).Our focus in the subsoil reconstruction, at these shallow depths, was on the geometry of sediments hosting aquifers and the impervious layers separating aquifers.The geological framework was reconstructed through several sections.Figure 2b shows an example of a section traversing the graben structure.
The lithostratigraphic correlations of wells and sections allowed synthesis of the conceptual geological/hydrogeological model illustrated in Figure 3.The following units were defined: Coastal Dune Unit (CDU): eolian and marine sand deposits.This extends from the coast to the western side of Pisa town [42].It contains an aquifer fed by meteoric and sea water.
Recent Alluvial Cover (RAC): clays and discontinuous sand bodies, with deposits from the Arno and Serchio rivers [40].These host small, discontinuous, phreatic, lenticular aquifers in sand.
Gravel Fans Unit (GFU): a succession of continental gravel deposits with interlayers of minor clay-rich sediments that stretch from the Monte Pisano fans to the Pisa plain [42].This highly permeable unit hosts several springs [33,43,44] and is a recharge area for the aquifers of the alluvial-sedimentary plain.The lithostratigraphic correlations of wells and sections allowed synthesis of the conceptual geological/hydrogeological model illustrated in Figure 3.The following units were defined: Coastal Dune Unit (CDU): eolian and marine sand deposits.This extends from the coast to the western side of Pisa town [42].It contains an aquifer fed by meteoric and sea water.
Recent Alluvial Cover (RAC): clays and discontinuous sand bodies, with deposits from the Arno and Serchio rivers [40].These host small, discontinuous, phreatic, lenticular aquifers in sand.
Gravel Fans Unit (GFU): a succession of continental gravel deposits with interlayers of minor clay-rich sediments that stretch from the Monte Pisano fans to the Pisa plain [42].This highly permeable unit hosts several springs [33,43,44] and is a recharge area for the aquifers of the alluvialsedimentary plain.
Sands and Gravels Unit (SGU): this unit consists of a succession of sands and gravels.These are Sands and Gravels Unit (SGU): this unit consists of a succession of sands and gravels.These are separated by a discontinuous Clay Unit (CU) that separates the SGU into the SGU-1 (upper) and SGU-2 (lower) aquifers.Strong lateral heterogeneities are observed in the SGU.This unit, with high permeability (hydraulic conductivity (K) equal to 10 −3 m/s in gravels and 10 −5 m/s in sands), hosts the multilayer aquifer system [33].
transitions into CDU sediments westward.Fine-grained sediments of the CSU underlie the SGU.The CSU covers the bedrock, which is heterogeneous in lithology and hosts thermal waters.To the east, the GFU (comprising coarse-grained deposits of Monte Pisano alluvial fans) is interdigitated with SGU sediments.Monte Pisano is a major meteoric recharge area for the confined multilayered aquifer hosted by the SGU (Figure 3).Within the area of interest, the pre-Neogene bedrock has been intersected by two drill holes: the Poggio 1 well, drilled for gas exploration, and the geothermal well San Cataldo 1 [17].The former, located about 7 km SSE of Pisa, encountered carbonate Mesozoic formations at 693 m below ground level (bgl), whereas the latter, drilled at the eastern periphery of Pisa township, intersected the same units at 570 m bgl.From geophysical data, the thickness of the carbonate formations has been estimated [17] to be around 1500 m.Other drill holes, Zannone and Tombolo 1, located just outside the area of interest, have also intersected the pre-Neogene basement (at 360 m and 2493 m bgl, respectively) encountering different formations, Scaglia and Macigno, respectively.The locations of the above drill holes are shown in Figure 4.
A detailed gravity survey, using roughly 1000 stations evenly distributed over an area of 220 km 2 , was used to reconstruct the geometry of the bedrock and to define the main buried structures affecting the Pisa sedimentary plain, as required for the subsequent geothermal modeling.All the gravity measurements were processed according to the standard procedures [45].The theoretical g values were computed using the GRS80 formula [46].The free-air correction used the formula that considers the station latitude and the second-order term of the height.The Bouguer correction, with the Bullard B term, reduces the infinite Bouguer slab to a spherical cap.Finally, the terrain was corrected, up to 30 km, with curvature correction beyond 10 km, using a right rectangular prism algorithm [47].The corrections were computed using a mass density of 2300 kg m -3 .
The resulting Bouguer anomaly map (Figure 4) is the summation of the gravity effects related to both deep and shallower anomalous bodies.Several strategies were introduced to separate the contributions of the deeper sources from the shallower ones.In particular, an upward continuation procedure, horizontal and vertical gradient analysis and least-squares residual anomaly separation, were applied to the Bouguer anomaly data.These units have the following spatial arrangement: the SGU is covered by RAC, which transitions into CDU sediments westward.Fine-grained sediments of the CSU underlie the SGU.The CSU covers the bedrock, which is heterogeneous in lithology and hosts thermal waters.To the east, the GFU (comprising coarse-grained deposits of Monte Pisano alluvial fans) is interdigitated with SGU sediments.Monte Pisano is a major meteoric recharge area for the confined multilayered aquifer hosted by the SGU (Figure 3).
Within the area of interest, the pre-Neogene bedrock has been intersected by two drill holes: the Poggio 1 well, drilled for gas exploration, and the geothermal well San Cataldo 1 [17].The former, located about 7 km SSE of Pisa, encountered carbonate Mesozoic formations at 693 m below ground level (bgl), whereas the latter, drilled at the eastern periphery of Pisa township, intersected the same units at 570 m bgl.From geophysical data, the thickness of the carbonate formations has been estimated [17] to be around 1500 m.Other drill holes, Zannone and Tombolo 1, located just outside the area of interest, have also intersected the pre-Neogene basement (at 360 m and 2493 m bgl, respectively) encountering different formations, Scaglia and Macigno, respectively.The locations of the above drill holes are shown in Figure 4.
A detailed gravity survey, using roughly 1000 stations evenly distributed over an area of 220 km 2 , was used to reconstruct the geometry of the bedrock and to define the main buried structures affecting the Pisa sedimentary plain, as required for the subsequent geothermal modeling.All the gravity measurements were processed according to the standard procedures [45].The theoretical g values were computed using the GRS80 formula [46].The free-air correction used the formula that considers the station latitude and the second-order term of the height.The Bouguer correction, with the Bullard B term, reduces the infinite Bouguer slab to a spherical cap.Finally, the terrain was corrected, up to 30 km, with curvature correction beyond 10 km, using a right rectangular prism algorithm [47].The corrections were computed using a mass density of 2300 kg•m −3 .
The resulting Bouguer anomaly map (Figure 4) is the summation of the gravity effects related to both deep and shallower anomalous bodies.Several strategies were introduced to separate the contributions of the deeper sources from the shallower ones.In particular, an upward continuation procedure, horizontal and vertical gradient analysis and least-squares residual anomaly separation, were applied to the Bouguer anomaly data.
Inversion was effected using a modified Bott's algorithm [48], which, besides other refinements, considers spatially variable hyperbolic depth-density functions, thus including the different compaction rates within the basin fill in the modeling process.Accounting for the lateral and depth variation of density resulted in very robust outcomes: e.g., without applying any constraints, the solution quickly converged to an upper surface for the bedrock that was within 2% of the San Cataldo 1 and Poggio 1 control points.Inversion was effected using a modified Bott's algorithm [48], which, besides other refinements, considers spatially variable hyperbolic depth-density functions, thus including the different compaction rates within the basin fill in the modeling process.Accounting for the lateral and depth variation of density resulted in very robust outcomes: e.g., without applying any constraints, the solution quickly converged to an upper surface for the bedrock that was within 2% of the San Cataldo 1 and Poggio 1 control points.
Two master faults dominate the deeper structure of the basin (Figure 5), both with downthrow on the western side.One is at the base of the Monte Pisano, almost rectilinear and parallel to the mountain range; the other is curved and more northerly oriented, and just west of the Pisa township.The sector within the two faults is a structural high, bordered on either side toward the south by significant lows.This setting could be interpreted as the result of block tilting, arising from the development of listric geometries along the easternmost fault.Two master faults dominate the deeper structure of the basin (Figure 5), both with downthrow on the western side.One is at the base of the Monte Pisano, almost rectilinear and parallel to the mountain range; the other is curved and more northerly oriented, and just west of the Pisa township.The sector within the two faults is a structural high, bordered on either side toward the south by significant lows.This setting could be interpreted as the result of block tilting, arising from the development of listric geometries along the easternmost fault.

3D Modeling
The geothermal GIS data were transferred to a 3D Petrel database in accordance with the framework defined in the subsoil conceptual model.A DTM (cell size 10 × 10 m), public data from the Tuscany Region website [49], was imported to elaborate the topographic surface.
As with the GIS, the Petrel project was set to the WGS_1984_UTM_Zone_32N (EPSG 32632) coordinate system, and 712 wells (Figure 2a) were selected and imported with their names, latitudes, longitudes, elevations, and depths.Down each borehole, profile well tops were located at the tops of lithological units, in accordance with the proposed conceptual model (Figure 3).Lithological sections were used as additional inputs.Well tops and additional inputs were processed through make/edit surface tools to create the surfaces and the isobaths of the different units.Surfaces were modeled using the "convergent interpolation algorithm" (Petrel ® , 2010, Schlumberger-Software Integrated Solutions, Parma, Italy), producing regular mesh grids of 50 × 50 m.Finally, each solid unit was developed from these surfaces.Together, these solid units form the 3D geological model of the Pisa basin (Figure 6).Elaborated by the inversion of gravity data.Coordinate system in meters.

3D Modeling
The geothermal GIS data were transferred to a 3D Petrel database in accordance with the framework defined in the subsoil conceptual model.A DTM (cell size 10 × 10 m), public data from the Tuscany Region website [49], was imported to elaborate the topographic surface.
As with the GIS, the Petrel project was set to the WGS_1984_UTM_Zone_32N (EPSG 32632) coordinate system, and 712 wells (Figure 2a) were selected and imported with their names, latitudes, longitudes, elevations, and depths.Down each borehole, profile well tops were located at the tops of lithological units, in accordance with the proposed conceptual model (Figure 3).Lithological sections were used as additional inputs.Well tops and additional inputs were processed through make/edit surface tools to create the surfaces and the isobaths of the different units.Surfaces were modeled using the "convergent interpolation algorithm" (Petrel ® , 2010, Schlumberger -Software Integrated Solutions, Parma, Italy), producing regular mesh grids of 50 × 50 m.Finally, each solid unit was developed from these surfaces.Together, these solid units form the 3D geological model of the Pisa basin (Figure 6).

3D Geothermal Modeling
The surfaces of the lithological units were imported as .txtfiles into PetraSim 5, the graphic interface for the TOUGH2 geothermal calculation code, to produce the geothermal model.

3D Modeling
The geothermal GIS data were transferred to a 3D Petrel database in accordance with the framework defined in the subsoil conceptual model.A DTM (cell size 10 × 10 m), public data from the Tuscany Region website [49], was imported to elaborate the topographic surface.
As with the GIS, the Petrel project was set to the WGS_1984_UTM_Zone_32N (EPSG 32632) coordinate system, and 712 wells (Figure 2a) were selected and imported with their names, latitudes, longitudes, elevations, and depths.Down each borehole, profile well tops were located at the tops of lithological units, in accordance with the proposed conceptual model (Figure 3).Lithological sections were used as additional inputs.Well tops and additional inputs were processed through make/edit surface tools to create the surfaces and the isobaths of the different units.Surfaces were modeled using the "convergent interpolation algorithm" (Petrel ® , 2010, Schlumberger -Software Integrated Solutions, Parma, Italy), producing regular mesh grids of 50 × 50 m.Finally, each solid unit was developed from these surfaces.Together, these solid units form the 3D geological model of the Pisa basin (Figure 6).

3D Geothermal Modeling
The surfaces of the lithological units were imported as .txtfiles into PetraSim 5, the graphic interface for the TOUGH2 geothermal calculation code, to produce the geothermal model.

3D Geothermal Modeling
The surfaces of the lithological units were imported as .txtfiles into PetraSim 5, the graphic interface for the TOUGH2 geothermal calculation code, to produce the geothermal model.Knowledge of bedrock geometry, obtained by gravity data inversion, also has a fundamental role in thermofluidodynamic simulation because bedrock acts as the heat source for the shallower portion of the investigated system.The simulation (Figure 7) therefore included the portion of the Pisa plain delimited by the two master faults and the areal limits of the gravimetric survey (Figure 5).
A grid with a horizontal cell size of 200 × 200 m was applied to the data.The vertical cell size varied by layer.Each unit, irrespective of its thickness, was subdivided into 10 cells using the top surfaces imported from Petrel (Figure 7).
Therefore, the grid was denser in the upper part of the model, where the reconstructed subsurface units had higher definition (Figure 7).This allowed detailed thermal characterization of the units and  S1, see Supplementary Materials).
Grid size selection was influenced by the maximum number of cells possible without imposing an excessive computational load (about 40,000 cells), considering also the iterations necessary to reach the best fit between the measured and calculated data.The grid could be modified locally to better resolve output temperatures in smaller selected areas.
Energies 2018, 11, x FOR PEER REVIEW 9 of 20 Knowledge of bedrock geometry, obtained by gravity data inversion, also has a fundamental role in thermofluidodynamic simulation because bedrock acts as the heat source for the shallower portion of the investigated system.The simulation (Figure 7) therefore included the portion of the Pisa plain delimited by the two master faults and the areal limits of the gravimetric survey (Figure 5).A grid with a horizontal cell size of 200 × 200 m was applied to the data.The vertical cell size varied by layer.Each unit, irrespective of its thickness, was subdivided into 10 cells using the top surfaces imported from Petrel (Figure 7).Therefore, the grid was denser in the upper part of the model, where the reconstructed subsurface units had higher definition (Figure 7).This allowed detailed thermal characterization of the units and aquifers (up to 200 m depth) most suitable for energy exploitation.Measured or literature-derived values for the density, permeability, porosity, thermal conductivity, and specific heat were associated with the cells in each lithostratigraphic unit (Table S1, see Supplementary Materials).
Grid size selection was influenced by the maximum number of cells possible without imposing an excessive computational load (about 40,000 cells), considering also the iterations necessary to reach the best fit between the measured and calculated data.The grid could be modified locally to better resolve output temperatures in smaller selected areas.
To run the simulation, input values of the temperature and pressure were assigned to each cell in the top and bottom layers of the model.Input pressure values followed a hydraulic pressure gradient.The average atmospheric temperature for the area is 16 °C (for the years 1942-2016, Pisa meteorological station) [50] and this was assigned as the constant input temperature at the top of the model (coinciding with the DTM surface) for the whole area (Figure 8).To run the simulation, input values of the temperature and pressure were assigned to each cell in the top and bottom layers of the model.Input pressure values followed a hydraulic pressure gradient.The average atmospheric temperature for the area is 16 • C (for the years 1942-2016, Pisa meteorological station) [50] and this was assigned as the constant input temperature at the top of the model (coinciding with the DTM surface) for the whole area (Figure 8).
Several tests were performed before assigning temperatures at the bottom of the simulation volume, i.e., the top of the sedimentary basin bedrock.This is critical because the bedrock temperatures control the heat transfer from the bedrock surface to the overlying sediments.Firstly, a constant temperature of 50 • C was assumed at the heating plate (the top surface of the bedrock) in accordance with temperatures measured inside the geothermal reservoir in the San Cataldo 1 well [17].Secondly, as the bedrock top surface has an irregular shape, and hence variable depth, temperatures were calculated based on the local geothermal gradient (50-60 • C•km −1 ) [16,17].According to calculation based on local geothermal gradient, temperatures at the bedrock top varies from about 32 • C in the North West portion of the simulation volume, to about 90 • C in the its South East portion (Figure 8).
After each test, the calculated temperatures were compared with direct temperature surveys in boreholes and with measured data from previous works [38].The results show agreement (deviation <6%) between the measured temperatures for aquifers and the simulated temperatures in the model cells at the corresponding depth in the aquifer.The best result was obtained using a local geothermal gradient value of 56 • C•km −1 , and this represents the steady-state condition of the model (Figure 9).Several tests were performed before assigning temperatures at the bottom of the simulation volume, i.e., the top of the sedimentary basin bedrock.This is critical because the bedrock temperatures control the heat transfer from the bedrock surface to the overlying sediments.Firstly, a constant temperature of 50 °C was assumed at the heating plate (the top surface of the bedrock) in accordance with temperatures measured inside the geothermal reservoir in the San Cataldo 1 well [17].Secondly, as the bedrock top surface has an irregular shape, and hence variable depth, temperatures were calculated based on the local geothermal gradient (50-60 °C •km -1 ) [16,17].According to calculation based on local geothermal gradient, temperatures at the bedrock top varies from about 32 °C in the North West portion of the simulation volume, to about 90 °C in the its South East portion (Figure 8).
After each test, the calculated temperatures were compared with direct temperature surveys in boreholes and with measured data from previous works [38].The results show agreement (deviation <6%) between the measured temperatures for aquifers and the simulated temperatures in the model cells at the corresponding depth in the aquifer.The best result was obtained using a local geothermal gradient value of 56 °C •km -1 , and this represents the steady-state condition of the model (Figure 9).Several tests were performed before assigning temperatures at the bottom of the simulation volume, i.e., the top of the sedimentary basin bedrock.This is critical because the bedrock temperatures control the heat transfer from the bedrock surface to the overlying sediments.Firstly, a constant temperature of 50 °C was assumed at the heating plate (the top surface of the bedrock) in accordance with temperatures measured inside the geothermal reservoir in the San Cataldo 1 well [17].Secondly, as the bedrock top surface has an irregular shape, and hence variable depth, temperatures were calculated based on the local geothermal gradient (50-60 °C •km -1 ) [16,17].According to calculation based on local geothermal gradient, temperatures at the bedrock top varies from about 32 °C in the North West portion of the simulation volume, to about 90 °C in the its South East portion (Figure 8).
After each test, the calculated temperatures were compared with direct temperature surveys in boreholes and with measured data from previous works [38].The results show agreement (deviation <6%) between the measured temperatures for aquifers and the simulated temperatures in the model cells at the corresponding depth in the aquifer.The best result was obtained using a local geothermal gradient value of 56 °C •km -1 , and this represents the steady-state condition of the model (Figure 9).The temperature outputs were then exported, in tabular format, with the coordinates of the cell centers.To facilitate data interpretation and visualization, these simulation results (more than 13,000 points with X, Y and Z spatial coordinates and temperature values) were imported into the Petrel software to generate a 3D subsoil temperature distribution model (Figure 10).IDW interpolation produced temperature maps at different depths.
The temperature outputs were then exported, in tabular format, with the coordinates of the cell centers.To facilitate data interpretation and visualization, these simulation results (more than 13,000 points with X, Y and Z spatial coordinates and temperature values) were imported into the Petrel software to generate a 3D subsoil temperature distribution model (Figure 10).IDW interpolation produced temperature maps at different depths.

Geothermal Applications
The underground thermal model allows the realization and visualization, using a 3D platform (Petrel) and GIS tools, of different temperature outputs integrating both measured and simulated data.As an example, Figure 11 is a 3D map of temperatures within the first confined aquifer of the SGU (at half thickness).Alternatively, slices at different fixed depths can be produced; Figure 12 shows temperature maps at depths of −50 m, −100 m, and −130 m a.s.l (i.e., 50, 100, and 130 m below sea level).The maps show that temperature inhomogeneity affects the subsoil of the Pisa sedimentary plain.Higher temperatures (Figure 12 5).The differences in subsoil temperatures may have several causes.Cold fresh water coming from the Pisa mountains meteoric recharge area (Figure 3) through the Monte Pisano alluvial fan aquifer reduces local temperatures in the confined multilayered aquifer through recharge of cold meteoric waters.Rising thermal waters in the northeast, near San Giuliano Terme, along master faults of the graben structure (Figure 5), could contribute to the higher temperatures at the eastern limit.The bedrock depth, which is shallower to the north (−500 m a.s.l) of the area and deeper in the southeast (from −1100-−1400 m a.s.l), and the consequent different thickness and type of overlying sediments also affects subsoil temperatures.
An "urban heat island" effect, already well described for shallow aquifers in densely inhabited towns [51], because of thermal pollution, may be responsible for the shallow (50 m a.s.l), positive temperature anomaly (Figure 12b) observed below Pisa downtown.

Geothermal Applications
The underground thermal model allows the realization and visualization, using a 3D platform (Petrel) and GIS tools, of different temperature outputs integrating both measured and simulated data.As an example, Figure 11 is a 3D map of temperatures within the first confined aquifer of the SGU (at half thickness).Alternatively, slices at different fixed depths can be produced; Figure 12   C at 130 m), corresponding with a structural low (Figure 5).The differences in subsoil temperatures may have several causes.Cold fresh water coming from the Pisa mountains meteoric recharge area (Figure 3) through the Monte Pisano alluvial fan aquifer reduces local temperatures in the confined multilayered aquifer through recharge of cold meteoric waters.Rising thermal waters in the northeast, near San Giuliano Terme, along master faults of the graben structure (Figure 5), could contribute to the higher temperatures at the eastern limit.The bedrock depth, which is shallower to the north (−500 m a.s.l) of the area and deeper in the southeast (from −1100-−1400 m a.s.l), and the consequent different thickness and type of overlying sediments also affects subsoil temperatures.
An "urban heat island" effect, already well described for shallow aquifers in densely inhabited towns [51], because of thermal pollution, may be responsible for the shallow (50 m a.s.l), positive temperature anomaly (Figure 12b) observed below Pisa downtown.
Simulation modeling provides a detailed thermal characterization of subsoil units and aquifers and so can support the choice of geothermal heat mining technology, allowing the shallow and low temperature geothermal resource of Pisa area to be utilized for heating and cooling through either open-or closed-loop GHPs.The resource present at higher depths, in permeable fractured carbonate formations in the areas of shallow bedrock, could be exploited for geothermal district heating (GeoDH) with direct use of thermal fluids and heat pumps.The thermal characterization of the subsoil (Figures 10-12), derived by the integration of measured and simulated temperatures, allows us to calculate the energy potential of the subsoil and hence the suitability of the different geothermal technologies.Simulation modeling provides a detailed thermal characterization of subsoil units and aquifers and so can support the choice of geothermal heat mining technology, allowing the shallow and low temperature geothermal resource of Pisa area to be utilized for heating and cooling through either open-or closed-loop GHPs.The resource present at higher depths, in permeable fractured carbonate formations in the areas of shallow bedrock, could be exploited for geothermal district heating (GeoDH) with direct use of thermal fluids and heat pumps.The thermal characterization of the subsoil (Figures 10-12), derived by the integration of measured and simulated temperatures, allows us to calculate the energy potential of the subsoil and hence the suitability of the different geothermal technologies.Simulation modeling provides a detailed thermal characterization of subsoil units and aquifers and so can support the choice of geothermal heat mining technology, allowing the shallow and low temperature geothermal resource of Pisa area to be utilized for heating and cooling through either open-or closed-loop GHPs.The resource present at higher depths, in permeable fractured carbonate formations in the areas of shallow bedrock, could be exploited for geothermal district heating (GeoDH) with direct use of thermal fluids and heat pumps.The thermal characterization of the subsoil (Figures 10-12), derived by the integration of measured and simulated temperatures, allows us to calculate the energy potential of the subsoil and hence the suitability of the different geothermal technologies.

Evaluation of the Geothermal Potential for Open-Loop GHP Systems
The heat stored in the phreatic or confined aquifers may be mined through doublets of wells (both production and reinjection wells) and a heat exchanger that represents the interface with the heat pump.This type of assemblage allows heat to be sustainably mined, without any loss or contamination of water, because all the water from the production well is reinjected restoring the hydraulic capacities of the aquifer.
The geothermal potential of the aquifers of the Pisa plain subsoil can then be estimated using the volumetric method equation [23,24]: in which Q (kJ) represents the potential heat content of the aquifer, while Qw and Qs (kJ) are the heat content stored in groundwater and sediments, respectively.V (m 3 ) is the aquifer volume, n is the porosity, Cw and Cs (kJ•m −3 •K −1 ) are the volumetric heat capacity of water and solids, respectively, and ∆T (K) is the difference between the undisturbed temperature of the aquifers/sediments and the temperature after heat mining.The volume of the SGU was obtained from the elaboration of the 3D model, using Petrel tools.The true volume of the lower part of the SGU, and hence its thermal potential, would be less than the modeled volume because it includes minor interlayered clay deposits.These were not modeled because of poor information in the deeper wells.Using this formula, the SGU (SGU-1 + SGU-2) aquifer system contains 3.81 × 10 13 kJ stored in groundwater and 5.76 × 10 13 kJ stored in solid rock, sands, and gravels (Table S2, see Supplementary Materials ).However, not all this energy can be effectively exploited, as open-loop GHP technology can use only the energy stored in groundwater.Assuming a recovery factor of 10% [52], the exploitable energy in the SGU aquifer is 3.81 × 10 12 kJ or 1.06 × 10 6 MWh.The ArcGIS tool "raster calculator" was used to apply the formula within each raster cell (50 × 50 m) and thus map the geothermal potential for open-loop systems down to a depth of 150 m (the depth of SGU-2) (Figure 13).To assess the value of the resource relative to potential demand, a heating period of 200 days, derived from the national subdivision of climatic zones, was used [53].Higher potentials (Figure 13) correspond to areas richer in sands, and consequently water, in the subsoil volume.Lower potentials occur in areas with more clay-rich lithologies.This helps in selecting the best areas for adoption of open-loop systems.
The heat stored in the phreatic or confined aquifers may be mined through doublets of wells (both production and reinjection wells) and a heat exchanger that represents the interface with the heat pump.This type of assemblage allows heat to be sustainably mined, without any loss or contamination of water, because all the water from the production well is reinjected restoring the hydraulic capacities of the aquifer.
The geothermal potential of the aquifers of the Pisa plain subsoil can then be estimated using the volumetric method equation [23,24]: in which Q (kJ) represents the potential heat content of the aquifer, while Qw and Qs (kJ) are the heat content stored in groundwater and sediments, respectively.V (m 3 ) is the aquifer volume, n is the porosity, Cw and Cs (kJ•m −3 •K −1 ) are the volumetric heat capacity of water and solids, respectively, and ΔT (K) is the difference between the undisturbed temperature of the aquifers/sediments and the temperature after heat mining.The volume of the SGU was obtained from the elaboration of the 3D model, using Petrel tools.The true volume of the lower part of the SGU, and hence its thermal potential, would be less than the modeled volume because it includes minor interlayered clay deposits.These were not modeled because of poor information in the deeper wells.Using this formula, the SGU (SGU-1 + SGU-2) aquifer system contains 3.81 × 10 13 kJ stored in groundwater and 5.76 × 10 13 kJ stored in solid rock, sands, and gravels (Table s2, see Supplementary Materials ).
However, not all this energy can be effectively exploited, as open-loop GHP technology can use only the energy stored in groundwater.Assuming a recovery factor of 10% [52], the exploitable energy in the SGU aquifer is 3.81 × 10 12 kJ or 1.06 × 10 6 MWh.The ArcGIS tool "raster calculator" was used to apply the formula within each raster cell (50 × 50 m) and thus map the geothermal potential for openloop systems down to a depth of 150 m (the depth of SGU-2) (Figure 13).To assess the value of the resource relative to potential demand, a heating period of 200 days, derived from the national subdivision of climatic zones, was used [53].Higher potentials (Figure 13) correspond to areas richer in sands, and consequently water, in the subsoil volume.Lower potentials occur in areas with more clay-rich lithologies.This helps in selecting the best areas for adoption of open-loop systems.

Evaluation of the Geothermal Potential for Closed-Loop GHP Systems
The heat stored in the subsoil can also be exploited by the means of vertical borehole heat exchangers (BHEs) connected to a heat pump; heat is extracted from the ground using closed-loop pipes in which a fluid (usually water) circulates.The vertical pipes are cemented with special high-conductivity grout to sediments or rocks at the wall of the borehole.The thermal conductivity, heat capacity and undisturbed temperature of the subsoil affect the thermal exchange between the ground and the BHEs.These parameters vary from place to place, and are linked to the geological, Energies 2018, 11, 1591 14 of 20 hydrogeological and geothermal settings of the site where the system is installed.An accurate knowledge of these parameters is fundamental for correct dimensioning of the system.The number and length of BHEs required for the closed-loop system depends on the thermal characteristics of the subsurface and on the required heating (or cooling) capacity.
Evaluation of the geothermal potential of closed-loop systems is commonly based on specific software (i.e., EED, Earth Energy Designer) or tabular data [54].In this work, the geothermal potential of closed-loop systems was evaluated applying the "G.POT" (Geothermal POTential) method [55].As above, we considered a heating period of 200 days [53].
Ground thermal properties were derived from the numerical simulation and the literature.The ground thermal conductivity (W•m −1 •K −1 ) was depth-averaged over the first 100 meters, which represents the average length of BHEs (Figure 14).
The heat stored in the subsoil can also be exploited by the means of vertical borehole heat exchangers (BHEs) connected to a heat pump; heat is extracted from the ground using closed-loop pipes in which a fluid (usually water) circulates.The vertical pipes are cemented with special highconductivity grout to sediments or rocks at the wall of the borehole.The thermal conductivity, heat capacity and undisturbed temperature of the subsoil affect the thermal exchange between the ground and the BHEs.These parameters vary from place to place, and are linked to the geological, hydrogeological and geothermal settings of the site where the system is installed.An accurate knowledge of these parameters is fundamental for correct dimensioning of the system.The number and length of BHEs required for the closed-loop system depends on the thermal characteristics of the subsurface and on the required heating (or cooling) capacity.
Evaluation of the geothermal potential of closed-loop systems is commonly based on specific software (i.e., EED, Earth Energy Designer) or tabular data [54].In this work, the geothermal potential of closed-loop systems was evaluated applying the "G.POT" (Geothermal POTential) method [55].As above, we considered a heating period of 200 days [53].
Ground thermal properties were derived from the numerical simulation and the literature.The ground thermal conductivity (W•m -1 •K -1 ) was depth-averaged over the first 100 meters, which represents the average length of BHEs (Figure 14).The mapping used the thermal conductivity values and thicknesses of alluvial deposits adopted for the numerical simulation.For the mostly clayey sediments a value of 1.8 W m -1 K -1 was used, while for the mostly sandy aquifer a value of 2 W•m -1 •K -1 was assigned.Horizontal variations in thermal conductivity are linked to different thicknesses of clayey and sandy sediments.The highest values of thermal conductivity occur in areas where clay sediments are subordinate to the sandy aquifer.The temperature distribution (°C) at 50 m depth (Figure 12b) derived from the numerical simulation was used as the undisturbed temperature of the ground.Assuming a BHE length of 100 m, the temperature at 50 m depth can be adopted as the mean temperature along the pipe.Volumerelated specific heat capacity values were derived from guidelines [54]: a value of 2.5 MJ•m -3 •K -1 was used.
A standard vertical BHE was used, according to prior methodology [55], for evaluation of the geothermal potential throughout the area.The parameters for this BHE, needed for the calculation of borehole thermal resistance and heat transfer with the ground, are listed in Table s3 (see Supplementary Materials).The mapping used the thermal conductivity values and thicknesses of alluvial deposits adopted for the numerical simulation.For the mostly clayey sediments a value of 1.8 W•m −1 •K −1 was used, while for the mostly sandy aquifer a value of 2 W•m −1 •K −1 was assigned.Horizontal variations in thermal conductivity are linked to different thicknesses of clayey and sandy sediments.The highest values of thermal conductivity occur in areas where clay sediments are subordinate to the sandy aquifer.The temperature distribution ( • C) at 50 m depth (Figure 12b) derived from the numerical simulation was used as the undisturbed temperature of the ground.Assuming a BHE length of 100 m, the temperature at 50 m depth can be adopted as the mean temperature along the pipe.Volume-related specific heat capacity values were derived from guidelines [54]: a value of 2.5 MJ•m −3 •K −1 was used.
A standard vertical BHE was used, according to prior methodology [55], for evaluation of the geothermal potential throughout the area.The parameters for this BHE, needed for the calculation of borehole thermal resistance and heat transfer with the ground, are listed in Table S3 (see Supplementary Materials).
Figure 15 shows the distribution of geothermal potential, for single standard vertical BHE, based on the use of closed-loop technology.Comparing the calculated kW values with the thermal conductivity (Figure 14) and undisturbed temperature of the ground (Figure 12b), it is immediately clear that the highest thermal potential (4.75 kW) occurs in the areas of highest thermal conductivity and temperature.The highest temperature values are indicated as being north of the Arno river.
The whole investigated area is very promising for the installation of closed-loop systems, since from a single vertical BHE, between 3.5 and 4.75 kW can be extracted for heating purposes.Generally, a GHP plant would include more than 1 vertical BHE, depending on the power of the GHP and heating requirements of the buildings.
Energies 2018, 11, x FOR PEER REVIEW 15 of 20 Figure 15 shows the distribution of geothermal potential, for single standard vertical BHE, based on the use of closed-loop technology.Comparing the calculated kW values with the thermal conductivity (Figure 14) and undisturbed temperature of the ground (Figure 12b), it is immediately clear that the highest thermal potential (4.75 kW) occurs in the areas of highest thermal conductivity and temperature.The highest temperature values are indicated as being north of the Arno river.The whole investigated area is very promising for the installation of closed-loop systems, since from a single vertical BHE, between 3.5 and 4.75 kW can be extracted for heating purposes.Generally, a GHP plant would include more than 1 vertical BHE, depending on the power of the GHP and heating requirements of the buildings.is the geothermal potential for a standard vertical BHE for heating only, a = 0.0701, T0 is the ground undisturbed temperature, Tlim is the minimum temperature of the heat carrier fluid, λ is the ground thermal conductivity, L is the length of the pipes, t'c = tc/365 where tc is the heating period, Rb is the borehole thermal resistance, u's = rb2/(4αts) and u'c = rb2/(4αtc) where rb is the borehole radius, α = λ/ρc is the thermal diffusivity of the ground and ts is the operating time of the system.For a more practice use, we derived kW values from the calculated MWh/y.

Evaluation of the Geothermal Potential for the Bedrock
As previously stated, the temperatures in bedrock range between 40° and 90 °C at depths from 500 m to 1400 m bgl.The bedrock is formed by carbonate formations, possibly hosting thermal fluids.These fluids may be utilized directly or increasing their temperature (before using fluids) through GHP or other technologies, in GeoDH plants.
The evaluation of this type of geothermal resource, deeper and hotter than those described above, was performed using the volumetric method [23,24] and the reconstructed bedrock unit.For economic reasons (well drilling cost), the calculation was limited to the area where carbonate formations hosting thermal waters are relatively shallow (<600 m, structural high of the bedrock.The volume, over an area of 53 km 2 , assuming a constant thickness of 200 m, is 10.6 km 3 .The parameters used for the bedrock geothermal potential calculation are summarized in Table s4 (see where QBHE (MWh/y) is the geothermal potential for a standard vertical BHE for heating only, a = 0.0701, T0 is the ground undisturbed temperature, Tlim is the minimum temperature of the heat carrier fluid, λ is the ground thermal conductivity, L is the length of the pipes, t'c = tc/365 where tc is the heating period, Rb is the borehole thermal resistance, u's = rb2/(4αts) and u'c = rb2/(4αtc) where rb is the borehole radius, α = λ/ρc is the thermal diffusivity of the ground and ts is the operating time of the system.For a more practice use, we derived kW values from the calculated MWh/y.

Evaluation of the Geothermal Potential for the Bedrock
As previously stated, the temperatures in bedrock range between 40 • and 90 • C at depths from 500 m to 1400 m bgl.The bedrock is formed by carbonate formations, possibly hosting thermal fluids.These fluids may be utilized directly or increasing their temperature (before using fluids) through GHP or other technologies, in GeoDH plants.
The evaluation of this type of geothermal resource, deeper and hotter than those described above, was performed using the volumetric method [23,24] and the reconstructed bedrock unit.For economic reasons (well drilling cost), the calculation was limited to the area where carbonate formations hosting thermal waters are relatively shallow (<600 m, structural high of the bedrock.The volume, over an area of 53 km 2 , assuming a constant thickness of 200 m, is 10.6 km 3 .The parameters used for the bedrock geothermal potential calculation are summarized in Table S4 (see Supplementary Materials).The geothermal energy potential contained within the first 200 m of the bedrock, limited to the area of gravimetric structural high, is equal to 8.30 × 10 14 kJ, or 2.31 × 10 8 MWh.Except for information from the S. Cataldo 1 and Poggio 1 deep wells, no data regarding the nature and permeability of the bedrock lithologies are available.Consequently, the calculated geothermal potential is to be considered only as approximate guidance.

Sustainability
The detailed 3D thermofluidodynamic modeling of subsoil defined by the steady-state conditions of the simulated volume (Figures 8 and 9), represents a snapshot of the Pisa plain system before any geothermal utilization.This is the undisturbed situation but would be affected by heat mining.To determine the level of sustainable use of the resource, it would be appropriate to undertake dynamic modeling, inserting in the 3D model and simulator shallow or deep doublets of wells or vertical BHE, and hence determine the impact on the subsoil resource of GHP and GeoDH exploitation.This would be an important result because it provides a means of regulating application of these technologies to avoid environmental impacts on the subsoil thermal state or on the water resource.

Conclusions
This paper demonstrates a methodology for assessment of renewable geothermal resources, based on the collection, storage, interpretation, and visualization of a large volume of data in a 3D environment, allowing optimal use of such information.
Application of this methodology to the alluvial-sedimentary geological environment of the Pisa plain, indicates that geological modelling and geothermal numerical simulations, can provide data on the amount of energy that shallow geothermal resources (sediments, rocks, and water) can furnish, including evaluation of the local shallow geothermal (low temperature) potential for GHP (open or closed loop systems) or for wider area heat mining from deeper and hotter geothermal aquifers.The study area was found to have high geothermal potential: the heat stored in water and rocks, up to 150 m depth, is about 2.5 × 10 7 MWh.Thematic maps provided spatially explicit estimation of the geothermal potential of different heat exchange technologies, supporting the optimal design of heating and cooling plants.
Low-temperature geothermal resources are ubiquitous in the shallow earth crust, making our approach (as exemplified in the case study) widely replicable.The approach produces information about local geothermal resources, guiding low temperature geothermal plant designers, geologists, and engineers, towards the most efficient solution for a specific site.In particular, knowledge of subsoil geology, derived from 3D reconstruction and a well-geodatabase, can facilitate the diffusion of GHP technology by supporting site selection and drilling of well doublets (for GHP open-loop systems) or vertical heat exchangers (for GHP closed-loop systems).
The methodology not only identifies favorable locations for a specific exchange technology, but also provides a robust conceptual model of the subsoil and supports the thermofluidodynamic modeling necessary for design of sustainable plants that preserve the thermal state of subsoil over the long-term.This approach can ensure environmental sustainability through best practice in the design and regulation of GHP and GeoDH systems, with their attendant environmental benefits.
Funding: This work was financially supported by "Geo4P" project (EnerGea, CoSviG and Regional Government of Tuscany) and University of Pisa, grant to A. Sbrana.
Clay and Sands Unit (CSU): clayey and sandy sediments below the SGU.Just a few deep boreholes intercept this unit [44].

20 Figure 5 .
Figure 5. 3D view of the bedrock upper surface (isobaths in m a.s.l) and master faults of the basin.Elaborated by the inversion of gravity data.Coordinate system in meters.

Figure 5 .
Figure 5. 3D view of the bedrock upper surface (isobaths in m a.s.l) and master faults of the basin.Elaborated by the inversion of gravity data.Coordinate system in meters.

Figure 5 .
Figure 5. 3D view of the bedrock upper surface (isobaths in m a.s.l) and master faults of the basin.Elaborated by the inversion of gravity data.Coordinate system in meters.
to 200 m depth) most suitable for energy exploitation.Measured or literature-derived values for the density, permeability, porosity, thermal conductivity, and specific heat were associated with the cells in each lithostratigraphic unit (Table

Figure 7 .
Figure 7. Section of the 3D model of the selected area.Underground vertical exaggeration: 10×.Topographic surface vertical 2×.RAC: recent alluvial cover; SGU: sand and gravel unit; CU: clay unit; CSU: clay and sands unit.

Figure 7 .
Figure 7. Section of the 3D model of the selected area.Underground vertical exaggeration: 10×.Topographic surface vertical 2×.RAC: recent alluvial cover; SGU: sand and gravel unit; CU: clay unit; CSU: clay and sands unit.

Figure 8 .
Figure 8. 3D visualization of input temperatures for top and bottom layers of the model: the upper layer coincides with the topographic surface (DTM), while the bottom layer is the top of bedrock.The simulation volume is the same as in Figure 7; see also Figure 5 for isobaths of bedrock top.

Figure 8 .
Figure 8. 3D visualization of input temperatures for top and bottom layers of the model: the upper layer coincides with the topographic surface (DTM), while the bottom layer is the top of bedrock.The simulation volume is the same as in Figure7; see also Figure5for isobaths of bedrock top.

Figure 8 .
Figure 8. 3D visualization of input temperatures for top and bottom layers of the model: the upper layer coincides with the topographic surface (DTM), while the bottom layer is the top of bedrock.The simulation volume is the same as in Figure 7; see also Figure 5 for isobaths of bedrock top.

Figure 9 .
Figure 9. Temperature simulation results in the steady state: (a) Color block plot; (b) N-S and E-W cross-sections.Model volume is the same as in Figures 7 and 8.
Figure 9. Temperature simulation results in the steady state: (a) Color block plot; (b) N-S and E-W cross-sections.Model volume is the same as in Figures 7 and 8.

Figure 10 .
Figure 10.3D point cloud of temperatures (from TOUGH2 simulation results) imported into the Petrel environment.Colored points indicate temperature values between the topographic (upper) and bedrock (lower) surfaces.Model volume is the same as in Figures 7 and 8. Coordinate system in meters.
) are observed, beneath Pisa city and north and south of the town, at depths of 50 m (23 °C), 100 m, and 130 m (27 °C), and correspond to the structural highs of the bedrock (Figures 4, 5 and 11).High temperatures (26-27 °C at 130 m of depth) are also observed at the eastern limit of the simulation area near the footwall of Monte Pisano.Lower temperatures occur in the central southeastern sector of the simulation area (18 °C at 50 m depth and 22 °C at 130 m), corresponding with a structural low (Figure

Figure 10 .
Figure 10.3D point cloud of temperatures (from TOUGH2 simulation results) imported into the Petrel environment.Colored points indicate temperature values between the topographic (upper) and bedrock (lower) surfaces.Model volume is the same as in Figures 7 and 8. Coordinate system in meters.
shows temperature maps at depths of −50 m, −100 m, and −130 m a.s.l (i.e., 50, 100, and 130 m below sea level).The maps show that temperature inhomogeneity affects the subsoil of the Pisa sedimentary plain.Higher temperatures (Figure 12) are observed, beneath Pisa city and north and south of the town, at depths of 50 m (23 • C), 100 m, and 130 m (27 • C), and correspond to the structural highs of the bedrock (Figures 4, 5 and 11).High temperatures (26-27 • C at 130 m of depth) are also observed at the eastern limit of the simulation area near the footwall of Monte Pisano.Lower temperatures occur in the central southeastern sector of the simulation area (18 • C at 50 m depth and 22

Figure 11 .
Figure 11.3D visualization of temperature of SGU-1 at aquifer mid-depth, in the Petrel workspace.Vertical exaggeration: 30×.Coordinate system in meters.

Figure 11 .
Figure 11.3D visualization of temperature of SGU-1 at aquifer mid-depth, in the Petrel workspace.Vertical exaggeration: 30×.Coordinate system in meters.

Figure 11 .
Figure 11.3D visualization of temperature of SGU-1 at aquifer mid-depth, in the Petrel workspace.Vertical exaggeration: 30×.Coordinate system in meters.

Figure 15 .
Figure 15.Geothermal potential (kW) for heating-only by geothermal closed-loop systems (for 1 standard vertical BHE).Depth: 100 m.The calculation follows [46] QBHE = [a (T0 − Tlim) λ L t'c]/[−0.619t'c log (u's) + (0.532 t'c − 0.962) log (u'c) − 0.455 t'c − 1.619 + πλRb]; where QBHE (MWh/y)is the geothermal potential for a standard vertical BHE for heating only, a = 0.0701, T0 is the ground undisturbed temperature, Tlim is the minimum temperature of the heat carrier fluid, λ is the ground thermal conductivity, L is the length of the pipes, t'c = tc/365 where tc is the heating period, Rb is the borehole thermal resistance, u's = rb2/(4αts) and u'c = rb2/(4αtc) where rb is the borehole radius, α = λ/ρc is the thermal diffusivity of the ground and ts is the operating time of the system.For a more practice use, we derived kW values from the calculated MWh/y.