Abstract
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., [1–3]). 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 ([4–6] 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 ([7–12] 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 [13–16]. 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 [14–18], (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.
Although alternative interpretations have been proposed regarding the origin of the thrust belt of the Apennines (see, for example, [19, 26–30]), 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 [31–33]. 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]).
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
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 [41–43]. 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., [10–85]), 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, 51–57] 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 [12–29, 31–58]. 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) [24–29, 31–59].
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).
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].
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]).
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 [24–29, 31–59]). This well reached a depth of 2909 m stopping in the middle of a thick pack of reflectors, and not penetrating the K horizon [16–24]. 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 [24–29, 31–64].
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 [9–11]. 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, [65–68, 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 (http://www.sinanet.isprambiente.it/it/sia-ispra/download-mais/dem20/view) and a geological map (i.e., http://www502.regione.toscana.it/geoscopio/geologia.html) 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 [9–11] 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., [11–18]), 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 (http://www.crop.cnr.it). 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.
(a)
(b)
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).
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 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.
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 [13–18].
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, 71–75].
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.
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).
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.
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 [21–29, 31–72]. 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 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.
(a)
(b)
(c)
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.
Acknowledgments
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 [65–68]. 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)