Table of Contents Author Guidelines Submit a Manuscript
Volume 2019, Article ID 8492453, 21 pages
Research Article

Synthetic Seismic Reflection Modelling in a Supercritical Geothermal System: An Image of the K-Horizon in the Larderello Field (Italy)

1Istituto per la Dinamica dei Processi Ambientali, CNR, Milano, Italy
2Istituto di Geologia Ambientale e Geoingegneria, CNR, Roma, Italy
3Istituto di Geoscienze e Georisorse, CNR, Firenze, Italy
4Istituto di Geoscienze e Georisorse, CNR, Pisa, Italy

Correspondence should be addressed to L. Petracchini; ti.rnc.gagi@inihccartep.oznerol

Received 27 July 2018; Accepted 4 November 2018; Published 28 March 2019

Guest Editor: Matteo Lupi

Copyright © 2019 R. de Franco 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.


We tested synthetic seismic reflection modelling along a seismic line (CROP-18A) in the geothermal field at Larderello (Italy). This seismic line is characterized by a discontinuous but locally very bright seismic marker, named K-horizon, which has been associated with various geological processes, including the presence of fluids at supercritical conditions. Geological and geophysical data were integrated in order to develop a 3D subsurface model of a portion of the Larderello field, where extremely high heat flow values have been recorded. In the study area, the K-horizon is particularly shallow and supercritical deep conditions were accessed at depth. The 2D model of the main geological units up to the K-horizon was extracted from the 3D model along the CROP-18A and used to generate the synthetic TWT stacked seismic sections which were then compared with the observed stacked CROP-18A seismic section. To build the synthetic sections, generated through the exploding reflector approach, a 2D velocity model was created assigning to each pixel of the model a constant P-wave velocity corresponding to the related geological unit. The geophysical parameters and the geological model reconstructions used in the modelling process derive from a multidisciplinary integration process including geological outcrop analogues, core samples, and geophysical and laboratory information. Two geophysical models were used to test the seismic response of the K-horizon, which is associated with (1) a lithological discontinuity or (2) a physically perturbed layer, represented by a randomized velocity distribution in a thin layer. For the latter geophysical model (i.e., the physically perturbed layer), we have tested three different scenarios changing the shape and the thickness of the modelled layer. Despite the reliable calibration implied by the use of homogeneous units, the seismic modelling clearly shows that the physically perturbed layer provides a better explanation of the reflectivity features associated with the K-horizon.

1. Introduction

Reducing uncertainties in geothermal exploration is essential in order to support the growth of this promising and sustainable energy source. This is particularly true when dealing with the exploration of supercritical geothermal systems in which the reservoir fluids are expected to be in a supercritical state (for pure water, and ) (e.g., [13]). These challenging unconventional high-temperature geothermal systems, usually associated with shallow magmatic intrusion, could provide significantly higher well productivities with respect to wells drilled in classic hydrothermal systems [3].

Exploration strategies of geothermal reservoirs may significantly benefit from synthetic seismic reflection profiles that enable prospective features on acquired seismic reflection data to be detected and the geological-geophysical interpretation and model reconstructions to be calibrated ([46] and references therein).

We studied the Lago Boracifero sector of the Larderello geothermal field in Italy, where superheated steam is currently exploited from both a shallow reservoir in sedimentary units and a deep reservoir in crystalline rocks ([712] and references therein). In the Larderello geothermal field, the presence of a deep and enigmatic seismic marker, the K-horizon, has long been identified [13] in both 2D and 3D seismic surveys. Rather than being represented by a single reflection event, the K-horizon is characterized by a strip of diffused reflectivity, which is spatially and vertically localized (about 100-500 ms TWT thick), associated with lateral variations in reflectivity and sometimes by a bright spot signature [1316]. The origin of the K-horizon has been interpreted as the manifestation of several processes including (a) the quartz α-β phase transition [17], (b) the brittle-ductile transition in extensional setting [1418], (c) a rheological variation associated with a lithospheric thrust [19], (d) natural hydrofracturing in the proximity of lithological or petrophysical changes [20, 21], (e) the development of thermometamorphic aureola at the top of a Quaternary granitic intrusion [10], and (f) presence of warm fluids at overpressure conditions. Following the results of the San Pompeo 2 well ( and extrapolated to the depth of 2930 m, [22]), which reached the vicinity of the K-horizon, the presence of supercritical fluids confined in a relatively thin layer has also been proposed to explain its high reflectivity (e.g., [10]). This hypothesis was used to explain teleseismic converted wave behaviour [23] and was recently explored by a deep drilling project [24].

Our study had four main aims: (1) the calibration of a 2D seismic lines CROP-18A, acquired in the Italian deep crust seismic project [25], through the use of a 3D geological-geophysical model, (2) the definition of a conceptual model of the area based on a rock physics model; (3) the modelling of the seismic response of the 3D geological-geophysical model with associated implications for processing the seismic reflection data, and (4) the definition of the seismic signature of reservoir rocks possibly hosting supercritical fluids.

2. Geological Setting

The Larderello field is located in the inner part of the Northern Apennines (Figure 1), a sector of the orogenic belt of the Apennines.

Figure 1: (a) Structural sketch map at a regional scale of the Northern Tyrrhenian Sea and Tuscany (modified after [77]). The main structural features and the related Pliocene-Quaternary basins are shown. The red polygon defines the study area. (b) Simplified geological map of the Larderello geothermal area (modified after [18]). The geographic coordinate system used is WGS84-UTM32N; EPSG code: 32632.

Although alternative interpretations have been proposed regarding the origin of the thrust belt of the Apennines (see, for example, [19, 2630]), the Apennines appear to have developed as the result of the westward-directed subduction of a limb of the Mediterranean Tethys along the retrobelt of the southwestward prolongation of the Alps. In this hypothesis, the westward-directed Apennines subduction gradually replaced the “east-ward” Alpine subduction [3133]. Following the fast eastward retreat of the west-directed subduction zone, the Meso-Cenozoic sedimentary covers situated along the Apulian-Adriatic passive margin were accreted within the Apennine orogen. Contractional deformations were then followed by back-arc extensional tectonics that cross-cut the thrust pile.

However, other authors have suggested a polyphase tectonic evolution of the Northern Apennines due to an eastward-migrating compressional tectonic at the front and a subsequent eastward-migrating extensional tectonic, which has affected the inner part of the orogenic belt since at least the Early Miocene ([34, 35]; Figure 2). A complex tectonic evolution during Miocene-Pleistocene with a prevalent contribution of compressive tectonics up to the Pleistocene Epoch has also been proposed (e.g., [36]).

Figure 2: Geological cross section of the study area (see the X-X′ trace in Figure 1 for the location) (modified after [10]). The contact metamorphic isograde, main seismic horizons (H and K), isotherms, and radiometric ages are reported. Note that this illustration represents one of the various geological frameworks proposed for the study area. Table 1 mentions the proposed alternative interpretations of the K-horizon.

The geophysical data available show that southern Tuscany is characterized by a shallow Moho discontinuity (20-25 km depth), a reduced lithosphere thickness due to the uprising asthenosphere [37], and an active extensional tectonic as inferred by borehole breakout analysis [38].

The tectonic pile characterizing the study area locally crops out along ridges separated by Neogene tectonic basins filled with continental and marine sediments from the Miocene to Quaternary in age (Figure 1).

The stratigraphic scheme adopted in this work (see also [10, 39, 40]) is made up, from top to bottom, of the following main units (Figure 3): (i)Neoautoctonous complex (Miocene-Quaternary)(ii)Ligurian complex (Jurassic-Eocene)(iii)Tuscan Nappe (Triassic to Miocene) and the Tectonic Wedge complex (Paleozoic-Triassic)(iv)Crystalline Basement (Pre-Cambrian?-Paleozoic; Pliocene-Quaternary?) formed by the Metamorphic Unit (composed of the Phyllitic complex, the Mica-schist complex, and the Gneiss complex) and the Intrusive complex

Figure 3: Conceptual sketch of the tectonic and stratigraphic setting of the Larderello area (modified after [10]). The surfaces modelled in the present work are the topographic surface, the Base of Neogene deposits, the top of the Tuscan Nappe plus the TWC, and the top of the Metamorphic unit. The surface representing the seismic K-horizon, not shown here, was extrapolated from the literature [11].

The Neoautoctonous complex is composed of marine, lacustrine, and continental deposits, characterized by conglomerates, sandstone, clays, marls, and evaporites. It is mainly localized in the tectonic basins, and it lies unconformably on a deformed substratum. The Ligurian complex is composed of the Jurassic ophiolite sequence of the oceanic crust and its Jurassic-Cretaceous sedimentary cover of the Ligurian units and by the Cretaceous-Oligocene turbidites of the Sub-Ligurian units [4143]. The Ligurian complex was thrust eastwards on the Tuscan Nappe during the Oligocene-Aquitanian. The Tuscan Nappe represents a sedimentary succession deposited from the Triassic to the Miocene on the paleomargin of the Adria Plate [43, 44]. The unit reflects the paleogeography evolution of the area passing from an evaporitic to shallow water marine environment followed by pelagic conditions. The lithostratigraphy of the Tuscan Nappe is characterized by evaporites, dolostone, and limestone (Triassic), followed by shallow water limestones (Rhaetic to early Jurassic), and by pelagic limestones, cherty limestones, and marls (Jurassic to Oligocene). During the Miocene, the marine deposits were covered by siliciclastic synorogenic deposits. In the Larderello area, the Tuscan Nappe is strongly delaminated, locally completely absent (e.g., [1085]), while a tectonic wedge complex (TWC) is present between the Tuscan Nappe and the Metamorphic unit. The TWC is formed by Paleozoic metamorphic rocks, Triassic metasiliciclastics and carbonates, and Upper Triassic evaporites of the Tuscan Nappe [10, 39, 40]. The Metamorphic Unit is composed of three main complexes [10]: (i) the Phyllitic complex made up above all mainly by metagreywacke (Cambrian-Devonian in age) and locally by carbonate–siliciclastic metasediments (Silurian-Devonian in age), (ii) the Mica-schist complex (Precambrian?-Early Paleozoic? in age), and (iii) the Gneiss complex (Precambrian?-Early Paleozoic?).

With regards to the intrusive complex, granites have been cored in several deep wells and range in age from 3.8 to 1.3 Ma [45]. Accordingly, some interpretations of the deep structure of the Larderello geothermal field propose the presence of very large batholites (e.g., [10, 11, 46]), and data suggest a still active magmatic emplacement with a partial melt occurring at depth of a few kilometers (e.g., [12] and references therein).

3. Geothermal Characterization

Larderello is a small town located in Tuscany (Italy) along the metalliferous hills, which is an important mining and industrial area (Figure 1). Geothermal electricity production began in the early 1900s. Superheated steams at Larderello and surrounding areas feed an installed capacity of 795 MWe [47]. The fields of Larderello, Travale, and Radicondoli belong to a unique deep system that covers an area of about 400 km2 [9].

Two reservoirs are used for power production. The shallow reservoir is hosted in sedimentary units and consists mostly of Mesozoic limestones and anhydrite dolostone. The deep reservoir is hosted in crystalline rocks, with temperatures exceeding 350°C [11, 48, 49] and consists of phyllite, mica-schist, gneiss, skarn, hornfelses, and granite.

Larderello is classified as a “convective intrusive geothermal play” [50], due to the presence of a plutonic heat source that feeds a wide hydrothermal system.

Geological and geophysical data and thermal numerical modelling support the hypothesis of still molten igneous intrusions ([12, 5157] and references therein), which represent the heat source of the geothermal system.

The study area is centred on the Lago Boracifero sector in the south-easternmost part of the field (Figure 1). Recently, research has been carried out to improve, through data integration and a multidisciplinary approach, the conceptual model of the Larderello geothermal area [1229, 3158]. New drilling technologies have also been tested to explore deep resources at supercritical (very high temperatures and pressures) conditions, by deepening an existing dry well (Venelle 2) [2429, 3159].

3.1. K-Horizon

2D and recent 3D seismic surveys have clearly detected two important seismic reflectors in the deep Larderello reservoir [13, 14, 16, 60]: (i) the H-horizon and (ii) the K-horizon (Figures 2 and 4).

Figure 4: A reflection seismic section from the literature (modified after [16]). This figure shows the two deep seismic markers, the H-horizon and K-horizon. The K-horizon upper boundary is shown in black.

The H-horizon is a discontinuous high-amplitude reflector and represents the current mining target, corresponding to highly productive intervals [10]. The K-horizon is a high-amplitude and locally bright-spot-type reflector, deeper and slightly more continuous than the H-horizon (e.g., [13, 54, 61]). In some areas, the K-horizon represents the upper boundary of a zone characterized by seismic reflectors with a lozenge-shape geometry [14, 15]. Although numerous wells have been drilled in the Larderello geothermal field (Figure 5 and Table S1 show wells used in the study area for this work), the K-horizon has never been reached and in the literature, its reconstruction derives mainly from seismic line interpretations [11, 16, 84, 86]. In 1982, the San Pompeo 2 geothermal well was drilled down to the proximity of the K-horizon, which in this area reaches its regional culmination at a depth of about 3-4 km (Figure 2). The sudden increase in temperature and unexpected blowout of the well prevented direct measurements at the bottom hole, but the temperature of >400°C and pressure of >24 MPa were extrapolated to the depth of 2930 m [22].

Figure 5: Simplified geological map of the study area (modified from The stratigraphic units have been merged to represent the geological units modelled in the present work. The wells used to constrain the geological model are shown. The dotted black line represents the segment of the CROP-18A seismic reflection line (Figure 6(b)) used in synthetic seismic reflection modelling (note in grey the continuation of the CROP-18A).

The origin and nature of the K-horizon are thus still highly debated in the literature and represent a key to a comprehensive understanding of the Larderello system. The models proposed in the literature to explain the nature of the K-horizon (Table 1) mostly imply the presence of fluids and a correlation with the isotherm (from [62]).

Table 1: Interpretation models proposed in the literature to explain the evidence of the K-horizon in seismic profiles in the Larderello geothermal system. The corresponding references are quoted.

Firstly, [14, 15, 62] interpreted the K-horizon as a brittle-ductile transition (BDT), relating the high-amplitude reflection to the fluids stored in the proximity of this interface. The authors included the activity of an extensional shear zone along the K-horizon in their interpretation, considering the reflector as a kinematically active rheological boundary. [18] also considered the role of listric shear zones as pathways for fluid migration from the K-horizon, thereby explaining the discontinuity of the reflector (local reduction of reflectivity).

Similar rheological variations with completely different tectonic implications were accounted in the interpretation of the K-horizon. For example, [19] related the horizon to a main lithospheric thrust.

A further contribution to the nature of the seismic reflector was made by [17], who also considered the volumetric thermal expansion of quartz, which is more abundant than the other minerals in the basement rock, during its α-β phase transitions. The K-horizon is presumed to be a layer of microfractured rocks roughly following the α-β quartz transition temperature which in this area is particularly shallow.

A different model was proposed by [10] who considered the K-horizon as the top of a young Quaternary granitic intrusion at the level of a thermometamorphic aureole hosting supercritical fluids.

Vanorio et al. and De Matteis et al. [20, 21] carried out a detailed analysis of various seismic parameters highlighting the occurrence of fluids at the level of the K-horizon which was attributed to natural hydrofracturing. In addition, the authors disregarded the role of BDT of the reflector and accounted for lithological or petrophysical changes immediately below the horizon.

Tinivella at al. [72] associated the high reflectivities of the K-horizon to warm fluid at overpressure conditions rising from below.

The well that has most closely approached the K-horizon is Venelle 2, which was redrilled and deepened in 2017 within the framework of the DESCRAMBLE project (see [2429, 3159]). This well reached a depth of 2909 m stopping in the middle of a thick pack of reflectors, and not penetrating the K horizon [1624]. However, the new data show that the accentuated seismic reflection at a depth of 2750-2800 m corresponds to a zone of increased thermal gradient (up to 0.3°C/m and temperature of about 507-517°C at 2900 m) and of decreased pressure fracturing. Pressure decrease is also associated with an increase in both the gas content in drilling fluids and the rate of penetration during drilling activities [63]. These results suggest that the temperature just below the K-horizon depth could be as high as 600°C, corresponding to the molten phase of granite, and that the observed increase in the thermal gradient may be the manifestation of a transient thermal state induced by the recent emplacement (<50 ka) of a granitic intrusion [2429, 3164].

4. Data and Methods

In order to generate the 3D geological-geophysical model of the study area, an integrated approach was used to combine subsurface and superficial data. A synoptic flowchart of the methodological workflow used is shown in Figure S1.

We used well data, published geological information, subsurface geological maps, and a seismic reflection profile [911]. The study area (Figure 5) is characterized by numerous deep geothermal wells, and data from 69 of them are available from public databases (Table S1, [6568, 87]). We carefully revised the original well head coordinates and converted them into our geographical system (i.e., WGS84-UTM32N; EPSG code: 32632). No information regarding their deviation surveys is available, hence all the wells were considered vertical. The stratigraphic data were revised according to the geological units described in Section 2.

The modelled geological units represent seismic units with specific seismic velocities (i.e., Neogene deposits, Ligurian complex, Tuscan Nappe, TWC, and Metamorphic Unit – Figure 3). The Tuscan Nappe and the TWC were treated as a single seismic unit as their seismic velocities are similar to the other geological units.

A digital elevation model (DEM) with a resolution of 20 m ( and a geological map (i.e., were used to better constrain the geological surfaces at a shallow level.

The geological map was projected onto the DEM surface and the geological boundary of the outcropping units. The bottom of the Neogene Unit and the top of the Tuscan Nappe plus TWC were digitized to define the emerging limits of the modelled surfaces. The area of the DEM characterized by outcrops of the Tuscan Nappe and TWC were selected and integrated into the input data as they define the emerging portion of the Tuscan Nappe plus TWC surface. In addition, geological subsurface maps [911] were georeferenced and digitized to constrain the Metamorphic Unit and the K-horizon at depth.

Although some faults have been interpreted in the area (e.g., [1118]), they were not modelled in the present work since their geometry is controversial and poorly defined by the available data. Therefore, the geological model represents continuous surfaces, clearly shaped and influenced by the paleogeographic setting and subsequent tectonic processes. However, even if this simplification may have an impact on the geological interpretation of the modelled surfaces, it has exiguous effects on our analysis concerning the seismic imaging of the overburden model and of the K-horizon in their geometry and geophysical parameters (e.g., [88]).

The dataset collected (well data, portion of DEM characterized by outcrops of the modelled units, geological boundaries, and digitized isobath maps) were then imported into the Petrel software (Schlumberger). Well data were integrated with superficial geological data to better constrain the trend of the modelled surfaces at shallow depths (Figure S2). For those geological units with few or no well data (i.e., Metamorphic complex and K-horizon surface), digitized isobath maps were used to better constrain the modelling process (Figure S3). The convergent interpolation method with a grid increment of 200 m in the - and -directions was applied to create the surface in Petrel.

The 3D geological-geophysical model was then used to run the synthetic seismic model. Modelled geological data were extrapolated from the 3D geological-geophysical model along a segment of the seismic line CROP-18A (from CDP 629 to CDP 941). The CROP-18A seismic reflection line (hereafter CROP-18A) is a seismic line acquired within the framework of the CROP project (Italian deep crust seismic project; [25]) available both in raw and processed stacked data versions, within the seismic database of the CROP consortium ( CROP-18A was acquired using dynamite sources spaced about 180 m apart and by means of a geophone cable with group spacing of 60 m (an asymmetric split-spread geometry with offsets −3780 m and 7620 m) and a gap of 150 m obtaining 3200% coverage. The sampling time and record length were 2 ms and 25 s, respectively. The stacked data section of the line (Figure 6(a)) was obtained using standard processing.

Figure 6: (a) CROP-18A stacked section. (b) Zoom of the CROP-18A stacked section of Figure 6(a). Stacked data in the CMP range from 629 to 941 were compared with the synthetic stacked sections. The trace of the studied segment of the CROP-18A section is reported in Figure 5.

The modelled geological surfaces were sampled along the CDPs of the CROP-18A. The collected points were subsequently interpolated with a bicubic spline, with 15 m of horizontal spacing.

To create the synthetic seismic section (Figure S1), a velocity model was imposed on the studied section by assigning to each pixel () the P velocities (constant velocity) adopted for the corresponding geological unit (Table 2) in the model calibration (see Section 5.3). A Gaussian velocity perturbation was adopted to explore an alternative hypothesis based on the physical rock-model for the K-horizon (see Section 5.4).

Table 2: Interval velocities applied for the modelled units and their ranges from literature data. Interval velocity of the geological unit below the K-horizon represents a major issue for the seismic exploration in the Larderello area, and it depends on the geological-geophysical interpretative hypotheses assigned to the K-horizon itself (see text for details).

The synthetic seismic stack sections of CROP-18A were generated using the exploding reflector seismoacoustic approach [69] which was developed in Matlab by the CREWES consortium [85] and was partly modified by us in this application. In detail, we interfaced the “afd_explode.m” script function of CREWES [85] with a main script which allowed the construction and/or the changes of the input velocity matrix deriving from the 3D geological-geophysical model. The wave modelling parameters were then set. Matlab scripts were developed to assess the difference between the synthetic outputs and the data and to manage the model and synthetic output renderings.

The exploding reflector approach provides a rapid calculation of the zero-offset synthetic stacked sections and helped us to calibrate and validate several geological and geophysical interpretative hypotheses regarding the geothermal reservoir model of the Larderello area.

With this approach, the zero-offset synthetic stacked section is obtained by locating the sources along all reflecting interfaces of the model and the receivers on the common mid points (CMPs). The exploding reflector section is almost equivalent to the zero-offset stacked section obtainable by standard CMP seismic processing in reproducing seismic wavefields not only by travel times but also by amplitudes.

In the present study, the positions of receivers were located at the CMP positions of the line between the CMPs 629 and 941, which were spaced 30 m apart, and a time window of up to 4 s of TWT was set (Figure 6(b)). The exploding reflector approach generates the seismograms for the P-wave velocity model obtained from the section studied.

The wavefield is propagated upward in depth using a finite difference algorithm and is then convolved with the input wavelet (Ricker wavelet with 25 Hz of central frequency) to produce the seismogram at the receiver.

The central frequency was selected through a spectral analysis on the raw data of the CROP-18A stacked section (Figure 7).

Figure 7: Raw data spectral analysis and central frequency estimate: mean spectrum on all traces (continuous line) and mean spectrum of selected trace with offsets greater than 2000 m (dashed line). Vertical lines indicate the central frequency value.

Figure 7 shows the results of two mean spectra calculated along the traces corresponding to shots of the studied segment of the CROP-18A (Figure 6). The first mean spectrum (continuous line in Figure 7) was obtained considering all the traces (i.e., 8064). For the second spectrum (dotted line in Figure 7), an offset selection (i.e., between 2000 m and 7770 m) was made to evaluate the signal spectral content for offsets characterized by clear deep reflected events in the shot gathers. In the calculation of both spectra, a time window of 5 s was used for each trace.

Both spectra showed a dominant peak at about 15 Hz. The estimated central frequencies for the two mean spectra were for the mean spectrum of all traces and for the mean spectrum of selected traces. These were calculated using , where () is the spectral point.

Although the two central frequencies and their related errors are quite similar, we chose 25 Hz as it is directly related to the signal propagated in deep structures.

The finite difference algorithm uses a nine-point approximation of the Laplacian operator and assumes the absorbing boundaries on the three sides of the model (bottom, right, and left). Time step modelling was set to 0.1 ms, and the Courant maximum number was ~0.2. The synthetic waveforms were sampled at 2 ms, similar to CROP-18A data.

Finally, our approach is essentially a kinematic calibration (fit of the 0-offset reflected arrival times) with the support of the quali-quantitative comparison of the complete wavefields.

To perform the calibration of the reconstructed 2D model along the CROP-18A, we then compared the main seismic horizons generated in the synthetic stacked sections processed for different models with the stacked data section along the line.

5. Results and Discussion

5.1. 3D Geological-Geophysical Model

The topography, the base of the Neogene deposits, the top of the Tuscan Nappe plus TWC, the top of the Metamorphic Unit, and the K-horizon were modelled (Figure 8). Through these surfaces, it has been possible to define the major geological units of the study area characterized by specific seismic velocities and hence potentially discernible along seismic sections.

Figure 8: Stratigraphic surfaces modelled with Petrel (Schlumberger). Well data and geological maps have been used to constrain the surfaces at depth.

The Neogene Unit and the Ligurian complex are not always present in the model (Figure 8). In particular, the Neogene Unit is not present in the south-eastern portion of the model which is characterized by outcrops of the Ligurian complex and the Tuscan Nappe plus the TWC (Figure 8). The Tuscan Nappe plus TWC is always present in the study area, and the top of this unit is well constrained by the wells available in the area. The Metamorphic Unit is located beneath the Tuscan Nappe plus TWC and at the bottom is confined by the K-horizon. The K-horizon was modelled on the basis of a published isobath map [11]. The modelled interfaces separating the main geological units were then extracted from the 3D geological-geophysical model along the studied segment of the seismic line CROP-18A (from CDP 629 to CDP 941) and subsequently converted into a time-domain geological-geophysical model.

5.2. Velocity Model

The geological units defined for the velocity model were the Neogene deposits, the Ligurian complex, the Tuscan Nappe plus TWC, the Metamorphic Unit, and the geological unit below the K-horizon. The velocity ranges defined for the aforementioned units derive from previous literature data and are reported in Table 2. Attributing a P-wave velocity for geological units may be difficult and generally has a range of values.

In order to obtain the seismic velocity model along the studied segment of the CROP-18A (from CMP 629 to CMP 941) a constant P-wave velocity (Vp) was assigned to each unit (Table 2). For the Neogene Unit and the Ligurian complex, we chose the arithmetic average of the interval velocities (i.e., 2700 m/s and 3850 m/s, respectively) commonly adopted in the literature [1318].

For the Tuscan Nappe plus the TWC, we adopted an interval velocity of 5700 m/s. This value, which is within the velocity ranges described in the literature (Table 2), was defined after a careful analysis of the data and taking into consideration the wide range of P-wave velocities reported in the literature (e.g., [13, 18, 76]) deriving from the considerable heterogeneity of these geological formations.

For the Metamorphic Unit, we adopted a velocity value of 5150 m/s, which is within the velocity range (i.e., 4400-5500 m/s) suggested by several authors (e.g., [6, 9, 13, 18, 76]). This velocity value was defined after analyzing the data observed in the Vp logs and in laboratory measurements of unfractured rocks from the Mica-schist complex and Gneiss complex in the depth range 1000-3500 m, as reported in the literature and summarized as follows [6, 7175].

Seismic reflectivity modelling of the deepest reflective horizons within the metamorphic units based on VSP measurement, logs, and AVO/AVA seismic analysis highlight that reservoir rocks are made up of fractured layers, with a variable thickness from one meter to tens of meters, showing a velocity variation of up to 15-30% with respect to the average velocity value of reservoir rocks [71, 72].

Vp logs of the Larderello-Travale area, at a depth of between 2400 to 3800 m, indicate a general negative asymmetric distribution around 20-30% of the Vp velocities with respect to a reference Vp value (5000-6000 m/s [73]). Furthermore, data logs indicate that the productive levels at the log scale are characterized by more complex velocity structures at short wavelengths (few meters) not resolvable with surface seismic. This means that the fractured layers are probably about 20 m thick with a velocity variation of up to about 5-10% [73]. The Vp laboratory measurements on core samples of the metamorphic unit (i.e., Mica-schist and the Gneiss complexes) of the Larderello area in a depth range of 1000-3800 m indicate a Vp anisotropy (differences in velocity in the vertical and horizontal directions of the core axis) of about 15% for pressures >40 MPa. In particular, for the core sample of the San Pompeo 2 well (2718.1 m), the Vp velocities in the vertical and horizontal directions are 5117 m/s and 5866 m/s, respectively [6, 74, 75]. Furthermore, the laboratory measurements [75] show that the linear-asymptotic trend, corresponding to fracture closure, in the pressure-velocity diagrams generally begins in the pressure range 100-150 MPa, and that Vp value variations of 15-20% and an anisotropy percentage of up to about 30% are observed in the pressure range 15-150 MPa. In addition, [74] report that both water-steam transition in the pore and high-temperature values produce only small changes in Vp measurements. They suggest that any major decreases in velocity could be associated with phase changes in fractures and fracture zones.

Due to its geological-geophysical complexity, the assignment of a P-wave velocity value for the K-horizon and for the unit below the K-horizon is even more complex than for the other units. In fact, it depends on the geological-geophysical hypotheses of the K-horizon itself. In order to reduce this complexity, there are two feasible hypotheses: (i) a sharp discontinuity which could be related for example to a lithological change or to an abrupt rheological transition and/or (ii) a perturbed layer with probable changes in the physical status of a portion of the geological unit and, therefore, the corresponding differentiation in the mechanical hydro-geophysical properties of rocks (i.e., confining, pore and effective pressures, temperature, porosity, and permeability) and type of inclusions and their petrophysical status (i.e., composition, salinity, and phases). This perturbed layer could be related, for example, to the presence of a thermometamorphic aureole, a mineral phase transition, or a highly fractured zone.

Previous studies, which focused on the processing and interpretation of the CROP-18A data and on the seismic velocity reconstruction from local earthquake tomography, suggest a Vp for the unit below the K-horizon ranging between 4300 m/s and 6400 m/s [6, 20, 21, 23, 61, 70, 72]. In addition to the uncertainties of each interpretation (about 10%), all these studies agree that the unit below the K-horizon should be assigned a high intrinsic variability of the Vp parameter, of up to about 15-40%, with respect to the Vp mean value adopted.

5.3. Model Calibration and K-Horizon Characterization by a Sharp Velocity Discontinuity

The geological-geophysical model along the CROP-18A seismic line was calibrated by superimposing the synthetic stacked sections on the seismic data stack for all the reflected events.

In order to simulate the reflection event corresponding to the K-horizon, a velocity contrast should be assigned to the unit below the K-horizon with respect to the upper Metamorphic unit. On the basis of the information discussed in Section 5.2 and on the hypotheses on the nature of the K-horizon and deep unit, we assigned a Vp of 4400 m/s and 5900 m/s to the unit below the K-horizon. These Vp values are close to the end members assumed for the units below the K-horizon [6, 20, 70, 72, 73] which correspond to a contrast of ±15% with respect to the Vp of the upper metamorphic unit.

In order to evaluate which of the two values should be used, we simulated the stacked sections for both models and velocities, and then we compared the seismic features (arrivals and polarities) of the reflection K-horizon with those observed in the data stacked section (8 and 9).

Figure 9 shows the starting velocity model (named model 1) obtained from the 3D geological-geophysical reconstruction. Here, the K-horizon interface was modelled with the smooth anticlinal shape proposed by [11]. Figure 10 shows the calibrated velocity model (named model 2) where the geometry of the K-horizon was modified to produce a better fit of synthetic stacked sections with respect to the section observed.

Figure 9: (a) Velocity model, named model 1-A, based on the 3D-geological reconstruction and imposing a velocity value of 4400 m/s for the geological unit below the K-horizon. (b) Velocity model, named model 1-B, based on the 3D-geological reconstruction and imposing a velocity value of 5900 m/s for the geological unit below the K-horizon. (a′) Modelled stacked sections for model 1-A (a) with P velocity of the unit below the K-horizon of 4400 m/s. (b′) Modelled stacked sections for model 1-B (b) with P velocity of the unit below the K-horizon of 5900 m/s. The modelled seismic response, plotted in positive variable area in red, is superimposed onto the processed stacked data, plotted in grey.
Figure 10: (a) Velocity model, named model 2-A, based on the 3D geological reconstruction obtained after the calibration with CROP-18A data and imposing a velocity value of 4400 m/s for the geological unit below the K-horizon. (b) Velocity model, named model 2-B, based on the 3D geological reconstruction after the calibration with CROP-18A data and imposing a velocity value of 5900 m/s for the geological unit below the K-horizon. (a′) Modelled stacked sections for model 2-A (a) with P velocity of the unit below the K-horizon of 4400 m/s. (b′) Modelled stacked sections for model 2-B (b) with P velocity of the unit below the K-horizon of 5900 m/s. The modelled seismic response, plotted in positive variable area in red, is superimposed onto the processed stacked data, plotted in grey.

Figures 9 and 10 report the velocity models (on top - a and b), the synthetic seismic responses (on the bottom - a′ and b′) for the two hypotheses on the P-wave velocity below the K-horizon, namely, 4400 m/s (a and a′ on the left) and 5900 m/s (b and b′, on the right). The synthetic response (in red variable area) is superimposed onto the stacked data (in grey wiggle and variable area) in order to compare the two zero-offset sections directly. This representation, in which the same plot parameters and trace normalization are used, enabled us to directly compare both data and synthetic responses in time domain and in relative normalized amplitudes.

The synthetic seismic responses indicate that model 2 is comparable with the seismic line drawing of all the main seismic reflection events detectable on the stacked CROP-18A line (Figure 10). For the K-horizon, the best seismic response of model 2 (Figure 10) was obtained by introducing a change in its original geometry as reconstructed in the 3D geological-geophysical model and reported in model 1 (Figure 9), consisting of a bulge-like structure located below the S. Pompeo 2 well area. In this area, the minimum depths of K-horizon for model 1 and model 2 are 3270 m and 2640 m, respectively.

The trace time lags obtained with trace cross-correlations between the data and synthetics signals calculated for the two velocities, in a window of 200 ms TWT including the K-event, have mean values of 5 ms and 18 ms and standard deviations of 40 ms and 73 ms, for 5900 m/s and 4400 m/s, respectively (Figure 10). The 5900 m/s simulation is characterized by the best fit with data and has a mean time lag comparable to the data sampling time of 4 ms. It reproduces the diffractions in the primary reflection between 1.2 s and 2 s TWT and the multiple events (Figure 10). Although this evidence would seem to indicate a positive contrast along the length of the K-horizon in our study area, the discontinuous pattern of the K-reflections means that there may be zones with a negative velocity contrast in other areas of Larderello.

It is interesting to note that the drilling operation for the S. Pompeo 2 well stopped at a depth of about 2767 m due to overpressure. This depth nearly corresponds to the top of the K-horizon, as modelled in our best fit of model 2.

In addition to the reliability of the reconstructed 3D geological-geophysical model, the results obtained by the seismic modelling indicates that (a)The seismic response for TWTs higher than 0.5 s is mainly influenced by the geometry of shallow units. This fact should be taken into account by including the effects of the Neogene deposits and/or the Ligurian complex in the static evaluation in order to optimize the seismic processing and the K-event focusing(b)Due to migration effects, the shape of the modelled K-horizon event in the distance range 2500-7500 m of the unmigrated seismic sections mainly depends on the geometry of the K-horizon in the model distance range of 4400-6700 m. The K-horizon in model 2 is characterized by a slightly more complex geometry than the reconstructed smoothed K-horizon of model 1 and introduces a more complex pattern of the K-related events into the time domain with unmigrated diffracted and multiple events also below the K-horizon. This pattern is in line with observations in the data stacked section. Consequently, when interpreting unmigrated K-events in the seismic section, there is a risk of assigning structural and physical meaning to nonexistent structures.

Figure 11 shows the time-migrated section using the calibrated model 2 and a Vp value of 5900 m/s below the K-horizon. The migrated section shows a clear improvement in the event focusing and a reduction in the diffracted events demonstrating the reliability of the calibrated model. In this context, we underline the high focusing of the edge-like structure between the CMPs 749 and 849 and the two events bordering the edge structure and dipping southward and northward. In particular, the northward dipping reflectors (right side of Figure 11) seem to continue upward to the surface at about CMP 700 and border the shallow basin structure located between the CMPs 700 and 820 (2500-4700 m of distance in the model, Figure 10). The migrated section shows, at a TWT greater than 2000 ms and below the edge structure, a diffuse reflectivity, and some events have an extension variable between several hundred meters to about 3 km (e.g., the event at 2250 ms to the right end of the studied segment).

Figure 11: (a) Velocity model, named model 2-B (see Figure 10(b)). (b) Studied segment of the CROP-18A stacked section (see Figure 6(b)). (c) The figure shows the time-migrated section of the CROP-18A data, reported in Figure 11(b), obtained through the calibrated velocity model shown in Figure 11(a). The time migrated section improves event focusing and reduces diffracted events. Note the high focusing of the bulge-like structure between CMPs 749 and 849 and the two events bordering the bulge structure and dipping southward and northward.
5.4. K-Horizon Alternative Hypotheses and Modelling: Physically Perturbed Layer

The results of our model calibration and line migration, and the processing and analyses of several seismic lines in the Tuscan geothermal area reported in the literature, highlight that the geological models and the related migration velocity models (above all for the deepest reflectors like H- or K-horizon) need to be more detailed in order to correctly migrate the reflectivity features observed in the area (i.e., [4, 6, 61]). In fact, the K-horizon exhibits some very particular reflective seismic features. It shows a lateral variation of reflectivity, and it is occasionally characterized by a bright spot signature. Also, instead of single reflection events, it shows a diffuse reflectivity spatially and vertically localized (generally extended 2-5 km in space, and 100-500 ms TWT in time) sometimes with a “lozenge” pattern below the top of the diffused reflective events [14, 15].

An alternative hypothesis is that the K-horizon is associated with a physically perturbed layer, as already claimed by [4]. This hypothesis is supported by the geological evidence derived from some outcrops on Elba, which are the remnants of an ancient and exhumed geothermal system considered representative of the deepest structural level of the Tuscan geothermal system, corresponding to the current level of the seismic K-horizon [77]. The results of geological studies carried out on the Elba outcrops suggest that the K-horizon is characterized by a permeability in the range 10−9 and 10−18 m2 and by fluid circulation. This fluid circulation has been inferred from the fluid inclusions of Elba, which reveal two main stages or groups of fluids: (a) hypersaline multiphase fluids (29-49 wt.% NaCl eq.) of lower temperature (<400°C) and (b) saline two-phase fluids (16-29 wt.% NaCl eq.) of a higher temperature (up to 600°C) [77].

Following this hypothesis, we modelled the K-horizon as a zone formed by a physically perturbed layer (PPL) characterized by a randomized P-wave velocity distribution. In our modelling hypothesis, the geological and physical properties inherent in the rock composition, density and shape of cracks, fluid characteristics, saturation conditions, and temperature and pressure conditions are reflected in the variations of velocity.

The high pore pressure in the region of a single fluid-filled pore helps to explain the resulting changes in the effective density and seismic velocities. Increased pressure within a pore tends to force the surrounding rock grains apart, thus tending to increase the pore volume. Effective density is a combination of the densities of the rock and fluid constituents of the porous medium. An increase in porosity results in a decrease in the effective density of the porous medium, as the rock is generally denser than fluid. A decrease in effective density tends to cause an increase in S-wave velocity. The effect of increased pore pressure on P-wave velocity can be seen by considering Wyllie’s time-average equation for P-wave velocity in porous, isotropic, fluid-saturated rocks under high pressure [78]. Mavko et al. [79] interpreted Wyllie’s equation as follows: the P-wave travel time through a fluid-saturated rock is equal to the sum of the travel time through the rock matrix and the fluid-filled pore.

In our model, we assume that the detection spatial scale of the velocity variations with seismic is comparable with the ¼ seismic wavelength (several ten of meters).

The pixel values () of the velocity perturbations inside the layer are assumed to be randomized with an asymmetric Gaussian velocity distributions according the velocity variations observed in data logs, laboratory measurements, and VSP and AVO-AVA analyses.

With an asymmetric distribution (Figure 12), the pixel velocities are set by fixing the highest Vp of the PPL equal to the one used in Section 5.3 for units beneath the K-horizon (i.e., 5900 m/s) and its negative variations are defined by the rock physical model described below.

Figure 12: The velocity distribution applied considering a physical perturbed layer (PPL) 100 m in thickness.

In order to insert a velocity perturbation that is consistent with a rock physical model, we assume that the confining pressure is equal to the lithostatic charge and that the pore pressure is equal to the hydrostatic charge, and the effective pressure is the difference between them. We assume a temperature of about 400°C and an effective pressure of about 30 MPa, which is in line with those derived from the measurement in the S. Pompeo 2 well [22]. A porosity of 3% and a density of 2700 kg/m3 were used. Density of 2700 kg/m3 is the average value for mica-schists according to [89], and very similar values have been used to model Larderello gravity data in [90]. The porosity of 3% is an average value for mica-schist as well, and a close confirmation was found in analogue mica-schist formation at Elba [77].

To calculate the negative velocity variations inside the PPL, we calculated the effective velocities using the scattering theory, considering a prolate (0.02 parameter) penny crack shape that characterizes the fluid inclusions [80]. To simulate the effect of the fluid in the pore and fracture space, we used the values for the velocity and density of brine (1000 m/s and 900 kg/m3) proposed by [81], as in other studies of the Larderello area by [2129, 3172]. Starting from a velocity of 5900 m/s for the unit below the K-horizon, we obtain the minimum velocity values of about 4500 m/s due to the presence of fluid-filled fractured layers, which represents the minimum value for the PPL velocity distribution (Figure 12). Furthermore, this value is about 11% less than the velocity value (5150 m/s) assigned to the overlying metamorphic rocks, causing negative velocity contrasts at the upper boundary of K-horizon. Finally, this value is in agreement with the values obtained by [72] and with laboratory velocity measurements in the same range of effective pressures (see Sections 5.2 and 5.3).

On the basis of the geological evidence observed on Elba analogue, three models were created (Figure 13): a model with the bulge structure below the San Pompeo 2 well characterized by a perturbed zone and two models with continuous perturbed layers thick 100 m and 500 m, respectively.

Figure 13: (a) Velocity model similar to the calibrated model in Figure 10 with a perturbed bulge. (b) Velocity model similar to the calibrated model in Figure 10 with a physical perturbed layer (PPL) of 100 m in thickness. (c) Velocity model similar to the calibrated model in Figure 10 with a physical perturbed layer (PPL) of 500 m in thickness. (a′) Modelled stacked section for the velocity model in (a). (b′) Modelled stacked section for the velocity model in (b). (c′) Modelled stacked section for the velocity model in (c). The modelled seismic response, plotted in positive variable area in red, is superimposed onto the processed stacked data, plotted in grey.

Figure 13 shows the synthetic stacked data obtained with the seismic responses (in red) superimposed on the stacked data (in grey).

The synthetic response for the PPL limited to the bulge structure (Figures 13(a) and (a′)) reproduces a K-horizon continuous event at the ends of the line (0-2600 m and 8000-9270), which derives from a positive velocity contrast. In this case, the K-horizon events are perturbed and lose continuity in the bulge area. On the other hand, there is a relatively continuous signal at the base of the bulge, characterized by a positive contrast. The response with a PPL of 100 m in thickness (Figures 13(b) and (b′)) is similar to the previous model although the K-horizon reflectivity at the extremities and at the base of the bulge becomes less continuous and reproduces a response that is more in agreement with the original stacked data. The response with a PPL of 500 m in thickness (Figures 13(c) and (c′)) shows a clear reflection at the base of PPL which is not clearly recognized in the stacked data. As expected, the higher the thickness of PPL, the higher the complexity of the reflection pattern below the K-horizon. The responses of model 2-B (i.e., the calibrated velocity model with a sharp discontinuity related, e.g., to a lithological change - Figure 10(b′)) and that with a PPL of 100 m in thickness (Figure 13(b′)) are quite similar, and they have the same best-fit of the primary reflection event of the K-horizon. The main difference is in slight perturbations of K-horizon event times, which are within the best-fit error limits, and in the lateral discontinuity of the signal amplitudes of the K-horizon event in the PPL model.

The amplitude lateral variation is the main reason for considering the PPL model as a more suitable physical model than model 2-B (Figure 10(b′)). The PPL may describe the amplitude features of the K-horizon reflection events observed in seismic explorations in the Tuscan geothermal region.

All the simulation results with the PPL (Figure 13) indicate a complex pattern of the multiple K-events (TWT greater than 2.3 s) with an emphasis on the lozenge feature.

We also investigated the influence of the central frequency of the Ricker source wavelet on the synthetic model response. Figure 14 reports the results considering the K-horizon considered as a PPL of 100 m in thickness (Figure 13(b)), in the TWT range 1-2.5 s, respectively, for 15 Hz, 25 Hz, and 40 Hz. These frequencies were selected on the basis of the results of the spectral analysis of raw data (Figure 7) which indicates the maximum peak at about 15 Hz, a central frequency of about 25 Hz, and its high frequency limit at about 40 Hz.

Figure 14: K-horizon synthetic response (100 m PPL model) for three different central frequencies of the source wavelet: (a) 15 Hz, (b) 25 Hz, and (c) 40 Hz.

The three responses indicate a different result passing from a relatively continuous reconstruction of the K-horizon for the 15 Hz source wavelet to a highly discontinuous reflector for the 40 Hz source wavelet. Although the high-frequency response is characterized by a general lateral high resolution, in reconstructing the PLL, the response generally has a low amplitude due to its sensitivity to low-velocity contrasts and interference.

6. Conclusions

The origin of strong amplitude seismic reflectivity concentrated in an almost continuous layer, named the K-horizon, below the vapor-dominated geothermal reservoirs of Larderello, Tuscany (Italy), has been debated for decades. In our paper, we contribute to the discussion by considering two hypotheses for the K-horizon: firstly, simulating the response of a sharp change in P-wave velocity at the depth of the K-horizon, and secondly, a physically perturbed layer (PPL) of variable thickness (representing the K-horizon) characterized by P-wave velocity perturbation.

A comparison of the two model reconstructions (Figures 10 and 13) clearly shows that the best fit is obtained by a PPL with a randomized P-wave velocity distribution with a thickness of 100 m. This is able to explain the reflectivity features and pattern of the seismic marker observed on several seismic profiles in the Larderello geothermal area.

The PPL modelled could be representative of a strongly fractured horizon filled by deep fluids possibly at supercritical conditions. The presence of a PPL is in agreement with the various hypotheses proposed in the literature for the K-horizon (Table 1), which may produce a P-wave velocity perturbation due to the change in the physical properties of the modelled rocks. Furthermore, the PPL hypothesis as proposed in this work is supported by the fact that the H-horizon, the current mining target in the study area made by fractured rocks with pressurized fluids, shows similar seismic features respect to the K-horizon. In this view, even the H-horizon could be characterized by a P-wave velocity perturbation due to the change in the physical properties of the modelled rock.

DESCRAMBLE project results ruled out the presence of fluids in the Venelle 2 well, but in contrast, the blowout of the San Pompeo 2 well clearly indicates the presence of deep overpressured fluids. This topic requires further research.

In addition to the possible presence of a PPL, this paper also highlights that (i)The deep reflection events are significantly influenced by the articulated shallow morphology of Neogene and Ligurian units. This influence should be taken into account in the seismic data processing of the line, and it poses the problem of the scale for the calculation of seismic statics. Data processing for the focusing of the deep reflecting horizons entails including the Neogene and the Ligurian units in the calculation of the statics(ii)The seismic modelling with a Vp value of 5900 m/s for the unit below the K-horizon is characterized by the best fit with data (Figures 10 and 13). These Vp values are in line with those reported in the literature for granites ([82, 83] and references therein; [84]). This may mean that there is a granitic intrusion below the K-horizon(iii)Different PPL thicknesses show a different seismic pattern of the reflectivity of the K-horizon. Hence, the lateral variation of reflectivity of the K-horizon, its bright spot signature, and the diffuse reflectivity could be a consequence of a lateral thickness variation in the PPL(iv)Our study also has implications regarding the processing and interpretation of seismic reflection data in geothermal areas(v)It confirms the difficulty in reconstructing suitable complex velocity models using only the velocity analysis from surface seismic and, as a consequence, in the migration of seismic lines. In this area, this issue is enhanced by the dual effect of static and the complexity of the K-horizon. The definition of the overburden heterogeneities (and anisotropy) using a high-resolution integrated seismic exploration approach up to 0.5 s TWT, with refraction tomography integrated with high-resolution reflection seismic, could limit the issue just to the complexity of deeper horizon response(vi)It suggests paying particular attention to the interpretation of seismic lines. In general, the response of the randomized model, due to constructive and destructive interferences, is not directly and univocally connected to the micro- and mesoscale structures of the reservoir. As a consequence, only the macroscale pattern can be recognized (see the reconstruction of the top of the physical transitional layer – Figure 13). Only the top of K-horizon can usually be detected, even though discontinuously. On the other hand, its internal structure may not be easily discernible due to the resolution decrease related to definition of a suitable local velocity model and to the central frequency of the signal (source and propagation)(vii)Reflectivity depends on the actual structure at the microscale, which cannot be directly reconstructed in detail from the processing of surface seismic data (vertical seismic profiling - VSP - has a better performance than surface seismic in this regard)(viii)An optimized seismic exploration strategy entails increasing the resolution of surface seismic, by integrating it with VSP data to visualize in detail the diffuse reflectivity zone detected by surface seismic

Data Availability

All the data used are public and available through the references and the website reported within the text.

Conflicts of Interest

The authors declares that they have no conflicts of interest.


We thank two anonymous reviewers, whose comments helped to improve and clarify this manuscript. This work was supported by the EU FP7-funded Integrated Methods for Advanced Geothermal Exploration (IMAGE) Project under grant agreement (608553).

Supplementary Materials

Figure S1: synoptic flowchart of the methodological workflow used. Figure S2: figure shows some of the data used to build the geological surfaces. In particular, well data of the Tuscan Nappe plus TWC and the area of the DEM characterized by outcrops of the Tuscan Nappe plus TWC are shown. Figure S3: figure shows the digitized subsurface map of the K-horizon imported in Petrel (Schlumberger). (a) View from above of the digitized data of the K-horizon. (b) Horizontal view of the digitized data of the K-horizon. Table S1: deep geothermal wells used to calibrate geological surfaces [6568]. The original well head coordinates have been converted into our geographical system (i.e., WGS84-UTM32N; EPSG code: 32632). No information regarding deviation survey is available. All the wells have been considered vertical. (Supplementary Materials)


  1. W. A. Elders, G. Ó. Friðleifsson, and A. Albertsson, “Drilling into magma and the implications of the Iceland Deep Drilling Project (IDDP) for high-temperature geothermal systems worldwide,” Geothermics, vol. 49, pp. 111–118, 2014. View at Publisher · View at Google Scholar · View at Scopus
  2. S. Scott, T. Driesner, and P. Weis, “Geologic controls on supercritical geothermal resources above magmatic intrusions,” Nature Communications, vol. 6, no. 1, article 7837, 2015. View at Publisher · View at Google Scholar · View at Scopus
  3. P. Dobson, H. Asanuma, E. Huenges, F. Poletto, T. Reinsch, and B. Sanjuan, Supercritical Geothermal Systems-a Review of Past Studies and Ongoing Research Activities, In 42nd Workshop on Geothermal Reservoir Engineering, 2017.
  4. M. Riedel, C. Dutsch, C. Alexandrakis, I. Dini, S. Ciuffi, and S. Buske, “Seismic depth imaging of a geothermal system in Southern Tuscany,” Geophysical Prospecting, vol. 63, no. 4, pp. 957–974, 2015. View at Publisher · View at Google Scholar · View at Scopus
  5. C. Schmelzbach, S. Greenhalgh, F. Reiser et al., “Advanced seismic processing/imaging techniques and their potential for geothermal exploration,” Interpretation, vol. 4, no. 4, pp. SR1–SR18, 2016. View at Publisher · View at Google Scholar · View at Scopus
  6. W. Rabbel, T. Jusri, D. Köhn et al., “Seismic velocity uncertainties and their effect on geothermal predictions: a case study,” Energy Procedia, vol. 125, pp. 283–290, 2017. View at Publisher · View at Google Scholar · View at Scopus
  7. A. Minissale, “The Larderello geothermal field: a review,” Earth-Science Reviews, vol. 31, no. 2, pp. 133–151, 1991. View at Publisher · View at Google Scholar · View at Scopus
  8. A. Barelli, G. Cappetti, and G. Stefani, “Results of deep drilling in the Larderello-Travale/Radicondoli geothermal area,” in Proceedings World Geothermal Congress 1995, pp. 1275–1278, IGA, Florence, Italy, 1995. View at Google Scholar
  9. R. Bertani, G. Bertini, G. Cappetti, A. Fiordelisi, and B. M. Marocco, “An update of the Larderello-Travale/Radicondoli deep geothermal system,” in Proceedings World Geothermal Congress, 2005, Antalya, Turkey, 2005.
  10. G. Bertini, M. Casini, G. Gianelli, and E. Pandeli, “Geological structure of a long-living geothermal system, Larderello, Italy,” Terra Nova, vol. 18, no. 3, pp. 163–169, 2006. View at Publisher · View at Google Scholar · View at Scopus
  11. P. Romagnoli, A. Arias, A. Barelli, M. Cei, and M. Casini, “An updated numerical model of the Larderello–Travale geothermal system, Italy,” Geothermics, vol. 39, no. 4, pp. 292–313, 2010. View at Publisher · View at Google Scholar · View at Scopus
  12. G. Gola, G. Bertini, M. Bonini et al., “Data integration and conceptual modelling of the Larderello geothermal area, Italy,” Energy Procedia, vol. 125, pp. 300–309, 2017. View at Publisher · View at Google Scholar · View at Scopus
  13. F. Batini, P. D. Burgassi, G. M. Cameli, R. Nicolich, and P. Squarci, “Combination to the study of the deep lithospheric profiles; “deep” reflecting horizons in Larderello-Travale geothermal field,” Memorie Società Geologica Italiana, vol. 19, pp. 477–484, 1978. View at Google Scholar
  14. G. M. Cameli, I. Dini, and D. Liotta, “Upper crustal structure of the Larderello geothermal field as a feature of post-collisional extensional tectonics (Southern Tuscany, Italy),” Tectonophysics, vol. 224, no. 4, pp. 413–423, 1993. View at Publisher · View at Google Scholar · View at Scopus
  15. G. M. Cameli, I. Dini, and D. Liotta, “Brittle/ductile boundary from seismic reflection lines of southern Tuscany (Northern Apennines, Italy),” Memorie Società Geologica Italiana, vol. 52, pp. 152–162, 1998. View at Google Scholar
  16. S. Ciuffi and M. Casini, “VSP and target definition,” in DESCRAMBLE project, Final Conference, Pisa, March 2018, View at Publisher · View at Google Scholar
  17. L. Marini and A. Manzella, “Possible seismic signature of the α-β quartz transition in the lithosphere of Southern Tuscany (Italy),” Journal of Volcanology and Geothermal Research, vol. 148, no. 1-2, pp. 81–97, 2005. View at Publisher · View at Google Scholar · View at Scopus
  18. A. Brogi, A. Lazzarotto, D. Liotta, and G. Ranalli, “Extensional shear zones as imaged by reflection seismic lines: the Larderello geothermal field (Central Italy),” Tectonophysics, vol. 363, no. 1-2, pp. 127–139, 2003. View at Publisher · View at Google Scholar · View at Scopus
  19. I. R. Finetti, M. Boccaletti, M. Bonini et al., “Crustal section based on CROP seismic data across the North Tyrrhenian–Northern Apennines–Adriatic Sea,” Tectonophysics, vol. 343, no. 3-4, pp. 135–163, 2001. View at Publisher · View at Google Scholar · View at Scopus
  20. T. Vanorio, R. De Matteis, A. Zollo, F. Batini, A. Fiordelisi, and B. Ciulli, “The deep structure of the Larderello-Travale geothermal field from 3D microearthquake traveltime tomography,” Geophysical Research Letters, vol. 31, no. 7, 2004. View at Publisher · View at Google Scholar · View at Scopus
  21. R. De Matteis, T. Vanorio, A. Zollo, S. Ciuffi, A. Fiordelisi, and E. Spinelli, “Three-dimensional tomography and rock properties of the Larderello-Travale geothermal area, Italy,” Physics of the Earth and Planetary Interiors, vol. 168, no. 1-2, pp. 37–48, 2008. View at Publisher · View at Google Scholar · View at Scopus
  22. F. Batini, G. Bertini, A. Bottai et al., “San Pompeo 2 deep well: a high temperature and high pressure geothermal system,” in Proceedings of the 3rd International Seminar on European geothermal update, vol. 3, pp. 341–353, FRG, Munich, Germany, 1983. View at Google Scholar
  23. N. Piana Agostinetti, A. Licciardi, D. Piccinini et al., “Discovering geothermal supercritical fluids: a new frontier for seismic exploration,” Scientific Reports, vol. 7, no. 1, p. 14592, 2017. View at Publisher · View at Google Scholar · View at Scopus
  24. R. Bertani, H. Büsing, S. Buske et al., “The First Results of the DESCRAMBLE Project,” pp. 12–14, Stanford University. View at Google Scholar
  25. D. Scrocca, C. Doglioni, F. Innocenti et al., Eds., “CROP Atlas: seismic reflection profiles of the Italian crust,” Memorie Descrittive della Carta Geologica d'Italia, vol. 62, pp. 1–194, 2003. View at Google Scholar
  26. E. Locardi and R. Nicolich, “Geodinamica del Tirreno e dell’Appennino centro-meridionale: la nuova carta della Moho. Memorie della,” Società Geologica Italiana, vol. 41, pp. 121–140, 1992. View at Google Scholar
  27. M. Boccaletti, G. Gianelli, and F. Sani, “Tectonic regime, granite emplacement and crustal structure in the inner zone of the Northern Apennines (Tuscany, Italy): a new hypothesis,” Tectonophysics, vol. 270, no. 1-2, pp. 127–143, 1997. View at Publisher · View at Google Scholar · View at Scopus
  28. C. Faccenna, F. Funiciello, D. Giardini, and P. Lucente, “Episodic back-arc extension during restricted mantle convection in the Central Mediterranean,” Earth and Planetary Science Letters, vol. 187, no. 1-2, pp. 105–116, 2001. View at Publisher · View at Google Scholar · View at Scopus
  29. G. Lavecchia, P. Boncio, N. Creati, and F. Brozzetti, “Some aspects of the Italian geology not fitting with a subduction scenario,” Journal of the Virtual Explorer, vol. 10, pp. 1–14, 2003. View at Google Scholar
  30. A. Brogi, A. Lazzarotto, D. Liotta, and G. Ranalli, “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. 1-2, pp. 60–80, 2005. View at Publisher · View at Google Scholar · View at Scopus
  31. C. Doglioni, E. Gueguen, P. Harabaglia, and F. Mongelli, “On the origin of west-directed subduction zones and applications to the western Mediterranean,” Geological Society, London, Special Publications, vol. 156, no. 1, pp. 541–561, 1999. View at Publisher · View at Google Scholar · View at Scopus
  32. E. Carminati, C. Doglioni, and D. Scrocca, “Alps vs. Apennines,” in 32° International Geological Congress, pp. 141–151, Italian Geological Society, Florence, Italy, 2004. View at Google Scholar
  33. G. Molli, “Northern Apennine–Corsica orogenic system: an updated overview,” Geological Society, London, Special Publications, vol. 298, no. 1, pp. 413–442, 2008. View at Publisher · View at Google Scholar · View at Scopus
  34. 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. 1-4, pp. 295–315, 1994. View at Publisher · View at Google Scholar · View at Scopus
  35. A. Brogi, “Neogene extension in the Northern Apennines (Italy): insights from the southern part of the Mt. Amiata geothermal area,” Geodinamica Acta, vol. 19, no. 1, pp. 33–50, 2012. View at Publisher · View at Google Scholar · View at Scopus
  36. M. Bonini and F. Sani, “Extension and compression in the Northern Apennines (Italy) hinterland: evidence from the late Miocene-Pliocene Siena-Radicofani Basin and relations with basement structures,” Tectonics, vol. 21, no. 3, pp. 1-1–1-32, 2002. View at Publisher · View at Google Scholar · View at Scopus
  37. G. Gianelli, “A comparative analysis of the geothermal fields of Larderello and Mt Amiata, Italy,” in Geothermal Energy Research Trends, H. I. Ueckermann, Ed., pp. 59–85, Nova Science Publishers, 2008. View at Google Scholar
  38. P. Montone, M. T. Mariucci, and S. Pierdominici, “The Italian present-day stress map,” Geophysical Journal International, vol. 189, no. 2, pp. 705–716, 2012. View at Publisher · View at Google Scholar · View at Scopus
  39. G. Gianelli, M. Puxeddu, and P. Squarci, “Structural setting of the Larderello-Travale geothermal region,” Memorie della Società Geoligica Italiana, vol. 19, pp. 469–476, 1978. View at Google Scholar
  40. E. Pandeli, G. Bertini, and P. Castellucci, “The tectonic wedges complex of the Larderello area (Southern Tuscany, Italy),” Bollettino della Societa Geologica Italiana, vol. 110, pp. 621–629, 1991. View at Google Scholar
  41. G. Nirta, E. Pandeli, G. Principi, G. Bertini, and N. Cipriani, “The Ligurian units of southern Tuscany,” Bollettino Società Geologica Italiana, vol. 3, pp. 29–54, 2005. View at Google Scholar
  42. M. Marroni and L. Pandolfi, “The architecture of an incipient oceanic basin: a tentative reconstruction of the Jurassic Liguria-Piemonte basin along the Northern Apennines–Alpine Corsica transect,” International Journal of Earth Sciences, vol. 96, no. 6, pp. 1059–1078, 2007. View at Publisher · View at Google Scholar · View at Scopus
  43. A. Brogi, D. Liotta, G. Ruggieri, E. Capezzuoli, M. Meccheri, and A. Dini, “An overview on the characteristics of geothermal carbonate reservoirs in southern Tuscany,” Italian Journal of Geosciences, vol. 135, no. 1, pp. 17–29, 2016. View at Publisher · View at Google Scholar · View at Scopus
  44. G. Bertini, F. M. Elter, and F. Talarico, “Evidenze di una fase estensionale pre-triassica nel complesso degli gneiss nell’area geotermica di Larderello (Toscana meridionale),” Studi Geologici Camerti, Volume Speciale, vol. 1, pp. 129–137, 1994. View at Google Scholar
  45. A. Dini, G. Gianelli, M. Puxeddu, and G. Ruggieri, “Origin and evolution of Pliocene–Pleistocene granites from the Larderello geothermal field (Tuscan Magmatic Province, Italy),” Lithos, vol. 81, no. 1-4, pp. 1–31, 2005. View at Publisher · View at Google Scholar · View at Scopus
  46. A. Arias, I. Dini, M. Casini, A. Fiordelisi, I. Perticone, and P. Dell’Aiuto, Geoscientific Feature Update of the Larderello–Travale Geothermal System (Italy) for a Regional Numerical Modeling, Proceedings world geothermal congress, 2010.
  47. R. Bertani, “Geothermal power generation in the world 2010–2014 update report,” Geothermics, vol. 60, pp. 31–43, 2016. View at Publisher · View at Google Scholar · View at Scopus
  48. F. Batini, A. Brogi, A. Lazzarotto, D. Liotta, and E. Pandeli, “Geological features of Larderello–Travale and Mt. Amiata geothermal areas (southern Tuscany, Italy),” Episodes, vol. 26, no. 3, pp. 239–244, 2003. View at Google Scholar
  49. G. Bertini, G. Gianelli, and A. Battaglia, Risultati ed interpretazione delle datazioni radiometriche (metodo 230Th/234U) dei campioni di minerali idrotermali presenti nelle rocce attraversate dai sondaggi geotermici (Larderello e Monteverdi) e negli affioramenti di rocce mineralizzate (Sassa e Canneto-Malentrata), ENEL-CNR-CISE internal report, Pisa, 1996.
  50. A. Santilano, A. Manzella, G. Gianelli et al., “Convective, intrusive geothermal plays: what about tectonics?” Geothermal Energy Science, vol. 3, no. 1, pp. 51–59, 2015. View at Publisher · View at Google Scholar
  51. J. E. Foley, M. N. Toksoz, and F. Batini, “Inversion of teleseismic travel time residuals for velocity structure in the Larderello Geothermal Field, Italy,” Geophysical Research Letters, vol. 19, no. 1, pp. 5–8, 1992. View at Publisher · View at Google Scholar · View at Scopus
  52. P. Baldi, S. Bellani, A. Ceccarelli, A. Fiordelisi, P. Squarci, and L. Taffi, “Correlazioni tra le anomalie termiche ed altri elementi geofisici e strutturali della Toscana meridionale,” Studi Geologici Camerti, Volume Speciale, vol. 1994-1, pp. 139–149, 1994. View at Google Scholar
  53. F. Batini, A. Fiordelisi, F. Graziano, and M. N. Toksöz, “Earthquake tomography in the Larderello geothermal area,” in Proceedings of the World Geothermal Congress, pp. 817–820, Florence, Italy, 1995.
  54. G. Gianelli, A. Manzella, and M. Puxeddu, “Crustal models of the geothermal areas of southern Tuscany (Italy),” Tectonophysics, vol. 281, no. 3-4, pp. 221–239, 1997. View at Publisher · View at Google Scholar · View at Scopus
  55. A. Manzella, G. Ruggieri, G. Gianelli, and M. Puxeddu, “Plutonic-geothermal systems of southern Tuscany: a review of the crustal models,” Memorie Società Geologica Italiana, vol. 52, pp. 283–294, 1998. View at Google Scholar
  56. F. Mongelli, F. Palumbo, M. Puxeddu, I. M. Villa, and G. Zito, “Interpretation of the geothermal anomaly of Larderello, Italy,” Memorie Società Geologica Italiana, vol. 52, pp. 305–318, 1998. View at Google Scholar
  57. A. Santilano, Deep Geothermal Exploration by Means of Electromagnetic Methods: New Insights from the Larderello Geothermal Field (Italy), [PhD Thesis], Politecnico di Torino, Italy, 2017.
  58. IMAGE Project website, “Integrated Methods for Advanced Geothermal Exploration,”
  59. DESCRAMBLE Project website, “Drilling in dEep, Super-CRitical AMBient of continenta L Europe H2020 project,”
  60. M. Casini, S. Ciuffi, A. Fiordelisi, A. Mazzotti, and E. Stucchi, “Results of a 3D seismic survey at the Travale (Italy) test site,” Geothermics, vol. 39, no. 1, pp. 4–12, 2010. View at Publisher · View at Google Scholar · View at Scopus
  61. F. Accaino, U. Tinivella, G. Rossi, and R. Nicolich, “Imaging of CROP-18 deep seismic crustal data,” Bollettino Società Geologica Italiana, vol. 3, pp. 195–204, 2005. View at Google Scholar
  62. D. Liotta and G. Ranalli, “Correlation between seismic reflectivity and rheology in extended lithosphere: southern Tuscany, Inner Northern Apennines, Italy,” Tectonophysics, vol. 315, no. 1-4, pp. 109–122, 1999. View at Publisher · View at Google Scholar · View at Scopus
  63. M. Cei and A. Fiorentini, “Resource characterization post drilling,” in DESCRAMBLE Project, Final Conference, Pisa, March 2018, http://www.descramble-h2020.Eu/index.Php/events/conferences. View at Publisher · View at Google Scholar · View at Scopus
  64. R. Bertani, “DESCRAMBLE project conclusive remarks,” in DESCRAMBLE Project, Final Conference, Pisa, March 2018,
  65. A. G. I. P. Mineraria, Temperature sotterranee, Inventario dei dati raccolti dall'AGIP durante la ricerca e la produzione di idrocarburi in Italia, San Donato Milanese (Milano), 1977.
  66. E.-A. ENEL and E. CNR, Inventario delle Risorse Geotermiche Nazionali—Indagine D’insieme Sul Territorio Nazionale, Ministero dell’Industria, del Commercio e dell’Artigianato (currently Ministero dello Sviluppo Economico), 1988.
  67. E. Barbier, S. Bellani, and F. Musmeci, “The Italian geothermal database,” in Proceedings World Geothermal Congress, pp. 3991–3995, Tohoku, Japan, 2000.
  68. “Ministry of Economic Development-DGS-UNMIG,”
  69. D. Loewenthal, L. Lu, R. Roberson, and J. Sherwood, “The wave equation applied to migration,” Geophysical Prospecting, vol. 24, no. 2, pp. 380–399, 1976. View at Publisher · View at Google Scholar · View at Scopus
  70. F. Batini, R. Bertani, B. Ciulli, A. Fiordelisi, P. Valenti, and E. G. Power, “Geophysical well logging: a contribution to the fractures characterization,” in Proceedings of the 27th Workshop on Geothermal Reservoir Engineering, pp. 28–30, Stanford University, Stanford, USA, 2002. View at Google Scholar
  71. G. M. Cameli, A. Ceccarelli, I. Dini, and A. Mazzotti, “Contribution of the seismic reflection method to the location of deep fractured levels in the geothermal fields of southern Tuscany (Central Italy),” in Proceedings World Geothermal Congress 2000, Kyushu-Tohoku, Japan, 2000. View at Publisher · View at Google Scholar · View at Scopus
  72. U. Tinivella, F. Accaino, G. Rossi, and R. Nicolich, “Petrophysical analysis of CROP-18 crustal seismic data,” Bollettino Società Geologica Italiana, vol. 3, pp. 205–211, 2005. View at Google Scholar
  73. M. Aleardi and A. Mazzotti, “A feasibility study on the expected seismic AVA signatures of deep fractured geothermal reservoirs in an intrusive basement,” Journal of Geophysics and Engineering, vol. 11, no. 6, p. 65008, 2014. View at Publisher · View at Google Scholar · View at Scopus
  74. B. Tarquini, M. Cameli, C. Donati, and R. J. Martin, “The effect of pressure, temperature, and lithology on seismic wave velocities for rocks from the Larderello geothermal field,” in Proceedings of the 1995 World Geothermics Congress, pp. 813–815, Florence, Italy, May 1995,
  75. G. P. Pinna, L. Burlini, D. Mainprice, and G. Gianelli, “Riflettività del basamento metamorfico di Larderello indicata da dati di petrofisica,” in 18° Convegno Nazionale GNGTS, Rome, Italy, 1999.
  76. F. Batini, S. Lisi, and E. Pinna, Contributo dei profili sismici in pozzo per l'individuazione di potenziali serbatoi geotermici all'interno del basamento metamorfico del campo di Bruciano (area di Larderello), 14° Convegno annuale del gruppo di geofisica della terra solida, Roma, 1995.
  77. M. Zucchi, A. Brogi, D. Liotta et al., “Permeability and hydraulic conductivity of faulted micaschist in the eastern Elba Island exhumed geothermal system (Tyrrhenian Sea, Italy): insights from Cala Stagnone,” Geothermics, vol. 70, pp. 125–145, 2017. View at Publisher · View at Google Scholar · View at Scopus
  78. M. R. J. Wyllie, A. R. Gregory, and L. W. Gardner, “Elastic wave velocities in heterogeneous and porous media,” Geophysics, vol. 21, no. 1, pp. 41–70, 1956. View at Publisher · View at Google Scholar
  79. G. Mavko, T. Mukerji, and J. Dvorkin, The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media, Cambridge university press, 2009. View at Publisher · View at Google Scholar
  80. J. G. Berryman, “Long‐wavelength propagation in composite elastic media I. Spherical inclusions,” The Journal of the Acoustical Society of America, vol. 68, no. 6, pp. 1809–1819, 1980. View at Publisher · View at Google Scholar · View at Scopus
  81. M. Batzle and Z. Wang, “Seismic properties of pore fluids,” Geophysics, vol. 57, no. 11, pp. 1396–1408, 1992. View at Publisher · View at Google Scholar · View at Scopus
  82. F. Birch, “The velocity of compressional waves in rocks to 10 kilobars: 1,” Journal of Geophysical Research, vol. 65, no. 4, pp. 1083–1102, 1960. View at Publisher · View at Google Scholar
  83. F. Birch, “The velocity of compressional waves in rocks to 10 kilobars: 2,” Journal of Geophysical Research, vol. 66, no. 7, pp. 2199–2224, 1961. View at Publisher · View at Google Scholar
  84. Y. Aizawa, K. Ito, and Y. Tatsumi, “Compressional wave velocity of granite and amphibolite up to melting temperatures at 1 GPa,” Tectonophysics, vol. 351, no. 3, pp. 255–261, 2002. View at Publisher · View at Google Scholar · View at Scopus
  85. CREWES-Consortium for Research in Elastic Wave Exploration Seismology, CREWES MATLAB Toolbox, University of Calgary,
  86. I. R. Finetti, “Basic regional crustal setting and superimposed local pluton-intrusion-related tectonics in the Larderello-M. Amiata geothermal province, from integrated CROP seismic data,” Bollettino Della Società Geologica Italiana, vol. 125, no. 1, pp. 117–146, 2006. View at Google Scholar
  87. E. Trumpy and A. Manzella, “Geothopica and the interactive analysis and visualization of the updated Italian National Geothermal Database,” International Journal of Applied Earth Observation and Geoinformation, vol. 54, pp. 28–37, 2017. View at Publisher · View at Google Scholar · View at Scopus
  88. C. Botter, N. Cardozo, S. Hardy, I. Lecomte, and A. Escalona, “From mechanical modeling to seismic imaging of faults: a synthetic workflow to study the impact of faults on seismic,” Marine and Petroleum Geology, vol. 57, pp. 187–207, 2014. View at Publisher · View at Google Scholar · View at Scopus
  89. A. S. Zappone and R. H. Bruijn, The Swiss Atlas of Physical Properties of Rocks (SAPHYR), Rapport Annuel, 2012.
  90. L. Orlando, “Interpretation of Tuscan gravity data,” Bollettino Società Geologica Italiana, vol. 3, pp. 179–186, 2005. View at Google Scholar