Simultaneous model-based inversion of pre-stack 3D seismic data targeting a deep geothermal reservoir, Northwest Hungary

The well-known, traditional way to extend P-wave acoustic impedance data between and beyond the well log locations is the post-stack inversion of seismic data usually available in the surroundings of the boreholes. A relatively new trend in the seismic exploration is based on the pre-stack inversion of the seismic CDP gathers providing both the P- and S-wave acoustic impedance sections (and volumes), as well as the estimated density data. This methodology is often called as simultaneous model-based inversion and can be utilized not only for hydrocarbon exploration, but it might also be a useful tool for the investigation of geothermal resources. In this study, we will compare the results of post-stack and pre-stack acoustic impedance inversions utilizing the same seismic volume. We will demonstrate and analyze the inverted attribute sections (and some of their derivatives) obtained by the pre-stack algorithm in detail. Finally, we will draw the conclusions about the lithological discrimination of the studied complex carbonate geothermal reservoir located in the pre-Cenozoic basement of the Little Hungarian Plain.


Introduction
In the late 1960s, the 'bright spot technology' resulted in several dry wells in the Gulf of Mexico. The reason was that high energy reflections can be caused not only by hydrocarbon content but also by other reasons, for instance magmatic or coal layers or even by the thin layer effect (Chopra and Castagna 2014). However, in the late 1970s, a more detailed analysis of the Direct Hydrocarbon Indicators (bright spots, flat spots, dim spots, polarity 1 3 change, high frequency attenuation) observed in the stacked seismic sections resulted more success in the exploration. The methodology became available after Taner et al. (1979) publication when they introduced the term of 'complex seismic trace'. Calculation of the reflection strength, instantaneous phase and polarity, and apparent frequency attributes based on the Hilbert-transform didn't need high-capacity computers; and they were calculated directly from the post-stack migrated seismic traces (without utilizing any well log data).
In the same year, another methodology was published by Lindseth (1979) to obtain P-wave acoustic impedance data from the reflected amplitudes of the post-stack seismic traces. This procedure became known as the Seislog technique and turned into a prevalent tool in the hydrocarbon exploration. This process utilized not only the post-stack seismic data, but also used well log information (in-situ sonic and density data) to build up initial model for the subsequent inversion process. This point can be regarded as the beginning of a new approach which was labelled as Seismic Lithology. The next step towards an even more reliable lithology (and porosity) estimation was Ostrander's (1984) fundamental publication. He proved that amplitude variations of the reflected signals with the source-receiver offset observed in the pre-stack CDP gathers were sensitive for the gas content of clastic sediments showing an anomalous behavior (increasing amplitudes with the offset). After introducing Amplitude Versus Offset (AVO) analysis into the hydrocarbon exploration, the number of productive wells increased significantly, and the method was introduced worldwide. There was a debate in the geophysical community whether AVO could work for fractured carbonates too (and not only for clastic sediments); but several authors proved that with successful case studies (such as Harvey 1993, Lynch et al. 1997, and Li et al. 2003. Finally, utilization of the AVO attributes (intercept, gradient, product, fluid factor, and Poisson's ratio change) became a reliable tool of the hydrocarbon indicators. Beyond the pre-stack seismic data, AVO inversion also requires well logs (sonic and density, and Full Wave Sonic at best). It needs very careful true amplitude seismic data processing before the lithological interpretation, and much higher computing power than the above-mentioned Hilbert and Seislog algorithms did. Goodway (2001) was the first who pointed out the importance of Lamé constants (incompressibility and rigidity) in lithology discrimination and  foreseen the possibility to estimate them from the pre-stack seismic amplitudes. To obtain the Lamé constants, we must perform a complex pre-stack workflow starting with AVO analysis, through the estimation of P-and S-wave acoustic impedances and the density and getting the incompressibility and rigidity finally. In this study, we intend to prove that the above workflow can be utilized not only for the hydrocarbon industry but also for the geothermal exploration. The reason is that the geological model (fluid saturated porous formations) and its rock-physical properties are very similar in both cases.
The first part of the above workflow (AVO analysis) was published by Masri and Takács (2022). In this paper, we will discuss our latest results around the possibilities of acoustic impedance inversion utilizing the same 3D seismic dataset and well logs.

Geology background and geophysical dataset
The study area (Fig. 1) is in the Hungarian part of the Danube Basin, which is the northwestern sub-basin of the Pannonian Basin (Sztanó et al. 2016). The sub-basin is characterized by vertical and lateral movements, pervasive extensional faulting, and high sandstone content in the Neogene sequence (Dolton 2006). Basement traps at the base Cenozoic include paleo-topography, truncations, and local porosity zones in the fractured and fissured traces of shear zones. In the western and northwestern part of the Danube Basin, the basement is dominated by high-grade metamorphic rocks. The Transdanubian Range, constituting the basement in the eastern and south-eastern part of the Danube Basin however, has only low-grade metamorphics and is dominated by various Mesosoic carbonate successions (Tari 1994;Fodor et al. 2003;Tari and Horváth 2010). The limit of the Transdanubian Range to the NW is the Rába Fault zone, a Miocene detachment fault system.
The geothermal gradient, and consequently the heat flow, is very high beneath the Danube Basin of the Pannonian Basin System (Békési et al. 2018). In the area of investigation, located on the north-west flank of the Transdanubian Range Unit, productive geothermal wells with higher effluent temperature than 100 °C can be expected only if they hit the basement (Cserkész-Nagy et al. 2020). The above described geological and geothermal settings provide an excellent environment for the exploration of hot water bearing reservoirs, mainly in the deeper parts, especially in the carbonate rocks of the pre-Cenozoic basement. Table 1 summarizes the geology and lithology revealed by a productive geothermal well located in the area of investigations.
A good selection of well logs was available in the above well, including measured sonic log which was essential to calculate acoustic impedance data at the location of the hole and to extend it along the properly processed pre-stack data of a 3D seismic cube.
In this study, we investigate the same 2 × 2 km portion of a larger 3D seismic dataset that was processed utilizing true amplitude data processing and AVO inversion to get reliable lithology information in our earlier paper (Masri and Takács 2022). Cropping the seismic volume was necessary to save running time for the time-consuming pre-stack Kirchhoff migration to get appropriate input data for the analyses. In the next chapters, we demonstrate two different kind of subsequent seismic acoustic impedance inversions focusing on the same geothermal reservoir located in the Triassic basement.

Acoustic impedance inversion methods using well log data
We will present the results of two different inversion processes (post-stack and pre-stack acoustic impedance inversions), as well as we will compare the obtained rock-physical attributes to assist the lithology and porosity discrimination in the study area. We will demonstrate the P-and S-wave acoustic impedances, density, and P-and S-wave velocity ratio attributes estimated based on the available 3D seismic data and well logs; and show the resulted parameters along the same vertical test profile.

Traditional post-stack acoustic impedance inversion
In Fig. 2, result of the traditional P-wave acoustic impedance inversion obtained by the post-stack method (Lindseth 1979) is shown assuming that the amplitudes of a stacked seismic trace are, on average, proportional to the reflection coefficients. Based on that approximation, Z(t) relative acoustic impedance variations can be calculated by the equation where Z 0 is the near surface acoustic impedance and R(t) is the reflectivity series down to the bottom of the investigated geological model. In other words, P-wave acoustic impedance can be estimated by integrating the deconvolved seismic trace, then exponentiation, and finally scaling by Z 0 (Lloyd and Margrave 2011). The procedure is applied on the stacked seismic traces providing 'high frequency' relative variations. The 'low frequency component' is calculated based on the product of sonic and density log values after applying a low pass filtering on the well log data. Finally, the 'absolute' seismic pseudoacoustic impedance trace is obtained by the sum of the relative variations and the low frequency component (Lindseth 1979), after scaling the result at the location of the well. This algorithm provides a single P-wave acoustic impedance (Z P ) volume. The calculation is fast and gives an appropriate result if we do not need S-wave impedance (Z S ) information and other useful rock-physical parameters -such as P-and S-wave wave velocity (V P and V S ) and density (ρ) for more detailed lithological investigations.

Simultaneous pre-stack acoustic impedance inversion
The significant difference between the above discussed post-stack acoustic impedance inversion and the pre-stack algorithm is demonstrated as follows. The pre-stack procedure provides not only a single P-wave impedance (Z P ) volume but also an S-wave impedance (Z S ) cube and other rock-physical parameters such as the P-and S-wave velocities (V P and V S ) and the density (ρ). After getting those characteristic lithology parameters, other very helpful indicators related to the Lamé parameters (λ-incompressibility and µ-rigidity) can be estimated by subsequent process called Lambda-Mu-Rho (LMR) transform (Goodway 2001, Russell and. Smith and Gidlow (1987) described the technology that can be used for converting the pre-stack reflected amplitudes of the CDP gathers into rock properties to estimate lithology and porosity. In practice, Amplitude Versus Offset (AVO) inversion can utilize several approximations to solve the complex Zoeppritz (1919) matrix equation (for example Aki Fig. 2 Inverted seismic P-wave impedance (Z P ) section calculated from post-stack amplitude data, and inserted P-wave impedance log at the location of the well where FR marks the top of the fractured zone and T1 marks its base based on the drilling reports. The magenta color denotes the highest acoustic impedance, and the green one denotes its lowest values and Richard 1980;Shuey 1985;Fatti et al. 1994). We note that our prior AVO inversion (Masri and Takács 2022) was performed using the Aki and Richard (1980) approach on the same pre-stack dataset.
In this study, Fatti et al.'s (1994) theory was utilized given as where R(θ) is the reflection coefficient versus incident angle (θ), Z P and Z S marks the P-and S-wave acoustic impedances, and V P and V S are the wave propagation velocities. The background V S /V P ratio should be known previously. Ma (2002) resolved this issue by replacing V S /V P with Z S /Z P . Other solutions based on the conjugate gradient method (Hampson et al. 2005;Russel et al. 2006) yield the P-wave and S-wave impedances, and the density parameters obtained from the pre-stack seismic data. This latter solution was applied in our present study (Hampson-Russell Software Services Ltd. 2007). Figure 3 shows the workflow of the simultaneous model-based inversion procedure based on the true amplitude pre-processed seismic gathers and the available well log data to obtain P-and S-impedance, and density parameters along the entire 3D seismic dataset. As we mentioned earlier, one of the great advantages of this procedure is that it can be followed by a subsequent algorithm to obtain even more indicative lithology and porosity indicators related to the Lamé parameters (incompressibility and rigidity).
The above inversion process is stabilized by the low frequency component of the available well log data applying low pass filtering. That low frequency information conditioned by the lithology is extended around the well locations using the interpreted horizons of the seismic sections (or volumes). The resulted start model is often called as 'low frequency or course layer' model. The 'high frequency' changes of the investigated parameters are calculated by the pre-stack inversion of the reflected amplitude of the CDP gathers. At the end, the final model is obtained with superimposing the high frequency variations to the low frequency (start) model after scaling the results at the location of the well logs. This procedure allows to extend the well-logging data (sonic and density) based on the pre-stack seismic reflections and provides the rock-physical parameters of Z P , Z S , and ρ in the surroundings of the well locations (Hampson-Russell Software Services Ltd. 2007). Certainly, and naturally, the resolution of the resulted seismic models is lower than the resolution of the well logs involved in the procedure.
(2) Fig. 3 Workflow of the simultaneous model-based pre-stack inversion procedure

Results and discussion
In this chapter we demonstrate the results obtained by the above discussed simultaneous pre-stack inversion and we analyze the calculated rock-physical attributes for mapping the hot water saturated fractured zone(s) of the study area. The first step of the procedure, similarly to any inversion algorithm, was the preparation of initial models (in our case V P , V S , and ρ) for the subsequent inversion process. The initial V P model was created by the spatial extrapolation of the low frequency (0-10 Hz) component of the measured sonic log utilizing a previously interpreted horizon as the top of reservoir around the 3D seismic volume. We note that the same start model was used for getting the result of the conventional post-stack inversion (see Fig. 2). For the pre-stack process, the initial S-wave velocity (V S ) and density (ρ) models were calculated from the low frequency P-wave velocity model (V P ) applying the Castagna et al. (1985) and Gardner et al. (1974) equations. The reason for those estimations was that measured Full Wave Sonic (V S ) and density (ρ) logs were not available in the study area at al. The lack of the measured Vs and ρ logs can affect a bit the accuracy of the final, inverted 3D rock physical models. However, the calculation of the final models was performed based on the measured pre-stack 3D seismic data (CDP gathers). The missing logs would have been used for only to get initial Vs and ρ models and for a calibration of the final models. If the initial models and the calibration were not exactly precise that does not affect the interpretation significantly. The reason is that mapping of the hot water saturated zone(s) was carried out based on the relative variations of the obtained rock physical parameters (decreasing or increasing). In this sense, the very precise absolute value of those models was not substantial.
The low frequency (0-10 Hz) initial P-wave acoustic impedance model with the inserted P-wave impedance log is presented in Fig. 4.
After setting up the low frequency initial models (0-10 Hz), high frequency variations (10-60 Hz) of the investigated parameters were calculated by the simultaneous inversion of the reflected amplitudes of the true amplitude processed CDP gathers. At the end, superposition of the initial models and the calculated relative changes provided the final models.
Before running the pre-stack simultaneous inversion on the whole 3D seismic dataset, an inversion analysis at the location of the well was helpful to verify the inversion parameters and optimize the scaling of the seismic amplitudes based on the well logs. Two tests were involved in this preliminary inversion analysis. The first one included linear regression to determine background relation between the rock properties to stabilize the inversion process. Cross-plots of the available and calculated well log data provided the following background equations for the study area: The second calculation of the inversion analysis performed at the well location was an estimation of a single scaler which was applied for the entire seismic dataset. This procedure sets up the RMS amplitude of the real seismic traces equal to the RMS amplitude of the synthetic seismic traces calculated from the well log data. Figure 5 shows the results: the measured and the inverted P-and S-wave impedance and the density logs; as well as the seismic synthetic, the real, and the error traces as a function of the incident angle. The wavelet to calculate the synthetic seismic traces was extracted from the real   Original (blue) and inverted (red) P-and S-wave impedance and density logs at the location of the well, as well as the seismic synthetic (red), real (black) and error (red) traces. The vertical axis is in twoway-time domain. The blue horizon in the seismic traces marks the top of reservoir predicted based on the conventional stack data and EOH means the end of hole seismic traces, using statistical wavelet estimations for the near and far offsets, in the whole depth range of the well logging. The correlation coefficient between the synthetic and real data was 0.654, which was not a too high value. However, the synthetic (red) and the real (black) traces showed very good similarity, and the error traces (red) didn't show any consistent flat events. Thus, we concluded that the inverted logs modeled the real data properly at the location of the well. We note that the depth range of the calculation window (1500-1850 m) included variable lithology from the Lower Pannonian, Miocene, and Triassic formations. Perhaps this was the reason for not getting any higher correlation value.
After the above local inversion analyses, the next step was to apply the simultaneous model-based inversion on the whole pre-stack 3D seismic dataset. Figure 6 demonstrates the inverted P-wave impedance model (Z P ) along a vertical section located closest to the productive well. Based on the color scale, it shows decreased Z P values inside the Triassic carbonates, between the top (FR) and the base (T1) of the known fractured zone. Note that low P-wave impedance value is one of the reliable indicators of the porous (fractured) rocks.
If we compare the above seismic attribute section obtained from pre-stack data ( Fig. 6) with the result of the preliminary post-stack inversion (Fig. 2), we can clearly see that the pre-stack algorithm provided very similar but more detailed image than the post-stack inversion did. Our recent finding agrees with Mallick's (2001) similar statement. However, neither of our pre-stack or post-stack P-wave impedance calculations was good enough to Fig. 6 Inverted P-wave impedance (Z P ) section calculated from the pre-stack seismic amplitude data and the available well logs. P-wave acoustic impedance log is inserted in the section at the location of the well distinguish the Triassic fracture zone (between the markers FR and T1) from the Triassic calcareous marl drilled beneath the T1 marker down to the end of hole (see Table 1). The reason could be that both the fractured dolomite and the calcareous marl can be characterized by low P-wave impedance values within the Triassic dolomite. We needed to create and study more rock-physical attributes to solve this issue. Figure 7 shows the inverted S-wave impedance (Z S ) section along the same test profile chosen very close to the well location. Looking at the inverted S-wave impedance data, it doesn't show any changes in the interval of the fractured zone (FR-T1). The reason is that the S-wave doesn't see the pores or fractures because it cannot travel in fluids, it can travel only in the rock matrix. Figure 8 presents the inverted density (ρ) section showing just a very moderate decrease in the fractured dolomite zone (between FR and T1). At this moment, we are not sure why the inverted density data doesn't detect the fracture zone more characteristically. One possible reason could be that there was no density log measured in the well, it was synthetized from the measured sonic log (Gardner et al. 1974). The other and very likely reason is that extracting density information from the seismic observations is very difficult because of its low sensitivity to the reflected amplitudes.
V P /V S velocity ratio is demonstrated along the same profile in Fig. 9. This is the most detailed image that we got in this study. It shows a zone of very low V P /V S values at the porous zone (between the FR and T1 markers), namely the dark green zone just below the top of the known reservoir (FR). Beneath the T1 marker, bit higher values appear with Fig. 7 Inverted S-wave impedance (Z S ) section calculated from the pre-stack seismic amplitude data and the available well logs. P-wave acoustic impedance log is inserted in the section at the location of the well yellow colors down to the bottom of the hole. Based on this V P /V S attribute section, we believe that inside the Triassic basement, the dark green colors indicate hot water bearing dolomite zones, and the yellow and red ones indicate marl and dolomite without fracturing (or only with a very low of porosity).
The explanation for the above-described lithology model, which fits very well to the drilling data (Table 1), is that the Poisson's ratio (and consequently the V P /V S ratio) values for marl are usually higher than they are for the fluid saturated dolomite as it can be seen in Fig. 10 (modified after Miles et al. 1989). We note that the Poisson's ratio and V P /V S ratio parameters are much more characteristic for the lithology than the P-wave velocity or the P-wave acoustic impedance.
Poisson's ratio (σ) is expressed by the following equation, where V P and V S are the Pand S-wave propagation velocities respectively: The V P /V S ratio attribute section presented in Fig. 9 also suggests significant spatial variations in the lithology of the Triassic basement. It is clearly seen, by the color scale, that the basement lithology is different at the two sides of the productive well. The reason for that can be the major fault zone, probably a strike slip zone, marked Fig. 8 Inverted density (ρ) section calculated from the pre-stack seismic amplitude data and the available well logs. P-wave acoustic impedance log is inserted in the section at the location of the well Fig. 9 Derived velocity ratio (V P /V S ) section calculated from the pre-stack seismic amplitude data and the available well logs. P-wave acoustic impedance log is inserted in the section at the location of the well Fig. 10 Poisson's ratio versus P-wave velocity for different lithology (redrawn after Miles et al. 1989) Note that Poisson' ratio is a nonlinear function of V P /V S and lower Poisson's ratio means lower V P /V S value. Color coding of the diagram is independent from the colors of the above seismic attribute sections with white segments. On the left side of the fault zone, we got moderate V P /V S ratio changes appearing with light green and yellow colors -very likely indicating friable dolomite and/or marl. On the right side of the fault zone, the variations are much more significant (dark green, yellow, and red colors) -detecting a cyclic series of high porosity dolomite, clayey dolomite with low porosity, and perhaps marl. Based on this image, we believe that the location of the productive well (drilled years before this study) was planned very carefully by focusing on the structural features of the prior (conventional) migrated stack volume. The recent study revealed a more detailed lithology model, and based on this, we can say that they drilled the edge of a complex geothermal reservoir of the Triassic basement at a highly permeable structural zone. That accurate well planning resulted in the very successful productive well supplying higher water temperature than 100 °C at the surface.
To accept or reject the lithology model discussed above, a subsequent Lambda-Mu-Rho (LMR) transformation will provide other estimated elastic parameters (incompressibility and rigidity). That procedure might give us even more helpful rock-physical information to verify the above model (Goodway 2001;Chopra and Castagna 2014).

Conclusion
P-and S-wave acoustic impedance, density, and V P /V S ratio attributes obtained by simultaneous model-based inversion of the available pre-stack seismic and well log data were very useful to investigate hot water bearing Triassic fracture zones and their surrounding lithology. The separation of the dolomitic fracture zones and marl formations was successful by the detailed images obtained around the area of investigation. Among the inverted rock-physical attributes, V P /V S ratio turned out to be the best for lithology discrimination and enabled to characterize the complex fabric of the pre-Cenozoic basement of the study area.
Our future task is to calculate other rock-physical attributes related to the Lamé parameters by the utilization of Lamda-Mu-Rho (LMR) transform. Those additional elastic parameters will contribute to a more detailed understanding of the lithology of the Triassic basement including the investigated complex geothermal reservoir.