Abstract

We have seismically characterized a Triassic-Jurassic deep geothermal sandstone reservoir north of Copenhagen, onshore Denmark. A suite of regional geophysical measurements, including prestack seismic data and well logs, was integrated with geologic information to obtain facies and reservoir property predictions in a Bayesian framework. The applied workflow combined a facies-dependent calibrated rock-physics model with a simultaneous amplitude-variation-with-offset seismic inversion. The results suggest that certain sandstone distributions are potential aquifers within the target interval, which appear reasonable based on the geologic properties. However, prediction accuracy suffers from a restricted data foundation and should, therefore, only be considered as an indicator of potential aquifers. Despite these issues, the results demonstrate new possibilities for future seismic reservoir characterization and rock-physics modeling for exploration purposes, derisking, and the exploitation of geothermal energy as a green and sustainable energy resource.

Introduction

In the future, geothermal resources may play a central role in supplying sustainable and green energy. For example, Danish geothermal resources in known sandstone aquifers can potentially sustain household heating needs for several centuries (Nielsen et al., 2004; Vosgerau et al., 2016). However, the ability to establish successful geothermal plants requires that the subsurface reservoir rock has a sufficient porosity and permeability. The application of appropriate exploration tools to constrain geologic risks and optimize site selection is, arguably, highly relevant when attempting to meet the needs of fragile geothermal budgets and uncertain profit models. Despite this, geothermal targets are often exclusively drilled based on conventional geophysical data interpretation techniques with little emphasis on more sophisticated quantitative reservoir characterization methods (Schmelzbach et al., 2016a).

Previous studies have demonstrated the advantages of geophysical exploration based on electrical resistance, electromagnetic, and gravity data for geothermal purposes (e.g., Nunziata and Rapolla, 1981; Pellerin et al., 1996; Barbier, 2002; Endo et al., 2018). However, for deep geothermal reservoirs, the seismic reflection method provides subsurface images with higher resolution. Previous geothermal studies based on seismic data contain several examples of structural and stratigraphic mapping of geothermal fields (Bujakowski et al., 2010; Casini et al., 2010; Lüschen et al., 2014, 2015). Furthermore, seismic attribute analysis allows the identification and characterization of high-permeability fractured zones (von Hartmann et al., 2012; Pussak et al., 2014; Khair et al., 2015), seismic tomography modeling (Muñoz et al., 2010), and vertical seismic profiling (Majer et al., 1988; Sausse et al., 2010; Place et al., 2011; Reiser et al., 2017). Several studies have applied passive microseismic measurements to monitor dynamic reservoirs with respect to microearthquakes induced by geothermal fluid stimulation (Dyer et al., 2008; Reshetnikov et al., 2015; Folesky et al., 2016; Schmelzbach et al., 2016b). Supplementary VP/VS ratio data from onshore multicomponent seismic measurements have also been used to discriminate between geothermal facies (Wei et al., 2014a, 2014b).

There are still, however, a scarce amount of studies that demonstrate the potential of combining modern seismic inversions with facies-dependent rock-physics modeling to improve our understanding of geothermal reservoirs in a quantitative manner. The hydrocarbon industry routinely uses such methods to seismically characterize reservoirs (e.g., Eidsvik et al., 2004; Buland et al., 2008; Grana and Della Rossa, 2010; Ba et al., 2013; Avseth et al., 2016), as well as for CO2 sequestration (e.g., Grude et al., 2013; Mallick and Adhikari, 2015; Grana et al., 2017). The benefits from such methods are equally valid for geothermal reservoirs because supplementary quantitative information about key parameters, such as porosity and permeability, can be used to improve site locations and predict geothermal reservoir performance. The objective of this study is, therefore, to extend seismic reservoir characterization to geothermal applications using a field example from an area north of Copenhagen in Denmark. We focus on predicting the facies and reservoir properties away from well locations based on seismic inversion data.

The data set available for this study was limited to a 2D seismic survey, one local well, and several regional wells that penetrate the same reservoir formations expected to occur at the prospect area. In addition, recent geologic and petrophysical analyses of the region provide valuable information that we use in the seismic inversion and rock-physics modeling (Nielsen et al., 2004; Mathiesen et al., 2010; Kristensen et al., 2016; Vosgerau et al., 2016, 2017; Weibel et al., 2017). For seismic reservoir characterization, we use previously reported methods (e.g., Avseth et al., 2005; Doyen, 2007; Dvorkin et al., 2014; Grana, 2018) that were implemented in a three-step workflow:

  1. 1)

    Rock-physics analysis: Well-log data were used to perform a feasibility study and calibrate facies-dependent rock-physics models.

  2. 2)

    Seismic inversion: Prestack seismic data were converted to absolute elastic parameters, such as the acoustic impedance and P-to-S velocity ratio, using a simultaneous amplitude-variation-with-offset (AVO) inversion.

  3. 3)

    Reservoir characterization: The results from steps 1 and 2 were combined to estimate facies and reservoir properties, such as porosity, permeability, and shale volume, in a Bayesian framework.

Uncertainties associated with intrinsic variabilities (e.g., the reservoir-property distributions beyond well control) and shortcomings in model theory or poor model calibrations were addressed using statistical methods. Probabilities and statistical estimators, which may provide valuable input for economic risk assessments and decision-making, represent the final facies and reservoir quality predictions (Avseth et al., 2005).

The primary geothermal target in this study is the sandstone-dominated upper Triassic-lower Jurassic Gassum Formation, which is currently being exploited for geothermal production and gas storage in Denmark (Røgen et al., 2015; Kristensen et al., 2016; Vosgerau et al., 2016, 2017). In addition, we evaluate the reservoir potential of a local lower Jurassic sandstone unit that overlies the Gassum Formation. Our results show that we are able to distinguish between several porous and clean sandstone bodies as potential aquifers and nonreservoir units within the target succession. Furthermore, our predictions indicate that a more promising reservoir potential exists within the lower Jurassic unit as compared with the Gassum Formation. Despite working with limited data coverage and quality, our predictions show a reasonably good match between the nearest well logs and geologic characteristics of the study area. Regional investigations prior to this study have also suggested that there are considerable lateral subsurface variations (Vosgerau et al., 2016, 2017), which emphasize the importance of supplementary quantitative seismic analysis.

In the following sections, we first introduce the key geologic settings and geothermal targets covered by our data set. We then describe the previously mentioned three-step workflow for seismic reservoir characterization. The reader should refer to the appendices for more theoretical details on the methodology and model specifications.

Data coverage and geothermal targets

Most of the Danish subsurface is characterized by deep sedimentary basins that contain extensive low-enthalpy geothermal resources (Balling et al., 2002; Mathiesen et al., 2010), which can potentially be exploited for district heating. Broad and uneven coverage of seismic reflection profiles and exploration wells, with varying data quality, constitute the primary subsurface data that are available for the entire country. This study focuses on geothermal exploration in an area near Hillerød in northeastern Zealand, Denmark (Figure 1). A 2D seismic survey performed in 2013 provides 3 km offset prestack, depth-migrated data and the nearest exploration well, Karlebo-1A, is located approximately 140 m from the nearest common-midpoint (CMP) gather along the eastern part of seismic line number 5 (see Figure 1). The seismic survey was designed for structural mapping and outlining potential geothermal reservoirs, whereas Karlebo-1A was drilled for hydrocarbon exploration.

The primary geothermal target in this study is the Gassum Formation, which is recognized for its excellent reservoir quality at several locations throughout Denmark. The sandstone-dominated lower Jurassic Reservoir Unit (LJRU) represents a secondary geothermal target, which lies above the Gassum Formation in the northeastern part of Zealand. The LJRU is a local phenomenon and forms the lower part of the regional cap rock, that is, the Fjerritslev Formation, which is dominated by marine mudstones and shales. The Gassum Formation and LJRU are similar with respect to lithology and are dominated by fine- to medium-grained, locally coarse, light-gray sandstones that alternate with mudstones, siltstones, and thin local coal seams (Michelsen et al., 2003). However, observations from the Karlebo-1A well (see Figure 1) indicate that the LJRU locally exhibits a slightly more homogeneous reservoir unit, whereas the Gassum Formation contains an increased number of interchanging thinly bedded shales. The reservoir temperature approximately ranges from 50°C to 65°C within the target interval (Poulsen et al., 2017).

Figure 2 shows a geologic interpretation of the geothermal targets and surrounding formations based on the seismic backdrop from line 5, with a superimposed Karlebo-1A projection. Obtaining seismic imaging of the geothermal targets is challenging, which is partly due to the overlying chalk package that generates interfering multiples and converted waves (Montazeri et al., 2018). The stratigraphic sequence interpretation of the composite succession, therefore, strongly guides the seismic interpretation. Certain mounded seismic characteristics within the geothermal targets can be observed and may represent possible seismic facies (Vosgerau et al., 2016), but we are unable to easily clarify whether these characteristics are related to variations in the geology or geophysical artifacts.

Rock-physics analysis

As a first step toward seismic reservoir characterization, we performed a feasibility study to verify whether the seismic properties are sensitive to the different facies and reservoir variations based on the regional well-log data (Waters and Kemper, 2014). With a positive outcome, we proceeded by calibrating the facies-dependent rock-physics models and used them to quantitatively describe the seismic inversion data in terms of the facies and reservoir properties. Bredesen et al. (2018) report such a regional feasibility study that we review and extend upon in this section. These analyses are primarily based on the Margretheholm-1A well, which we assume represents a close reservoir analog to the prospect area and is the only adjacent well available with a complete suite of logs.

To obtain an initial understanding of the seismic response of the target geothermal reservoirs, Figure 3 plots the acoustic impedance (Ip) against the P-to-S velocity ratio (VP/VS) data in a crossplot. Colorized dots and black “x” symbols represent clean (shale volume <30%) and shaly (shale volume >30%) sandstones from the LJRU and Gassum Formation and the cap rock (i.e., the Fjerritslev Formation), respectively. Kernel-density-estimated probability density functions (PDFs) (Silverman, 1986; Grana, 2018) were superimposed onto the data and projected onto each axis to form marginal distributions that describe the corresponding data histograms. Whereas the clean and shaly sandstones exhibit overlapping and nonunique properties, the two sandstone facies can be segregated from the cap rock due to larger Ip and lower Vp/VS values. To obtain a reference facies classification profile, we used linear discriminant analysis based on porosity and shale volume logs (Grana et al., 2017) to segregate (1) clean sandstone with high reservoir quality, (2) shaly sandstone with intermediate-to-low reservoir quality, and (3) nonreservoir facies.

Ultimately, we are interested in facies discrimination based on seismically inverted parameters, such as Ip and VP/VS, at a distance from well locations. Therefore, for the feasibility test, we used a Bayesian facies classification approach (Avseth et al., 2005; Doyen, 2007; Grana, 2018) to verify if the Ip and VP/VS parameters are sensitive to variations in the facies. The reservoir net-to-gross volumes were used as a priori information, with a 30%, 20%, and 50% probability for the clean sandstone, shaly sandstone, and nonreservoir, respectively. Figure 4a4e shows a suite of well logs from Margretheholm-1A, whereas Figure 4f4i shows the corresponding profiles for the reference and Bayesian-classified (or predicted) facies. To approach the seismic resolution, a sequential Backus averaging was used in a running window corresponding to one-fourth of the seismic wavelength to upscale the log data (the red stippled curves) and facies profiles (Avseth et al., 2005). The seismic wavelength was estimated to be 35 m using a center frequency of 50 Hz and an interval velocity of 3500 m/s extracted from the stacking velocities around the reservoir interval. In general, the upscaled profiles favored the shaly sandstone facies due to effective smoothing of the petrophysical and seismic properties. Nevertheless, the predicted and reference facies profiles match quite well at both scales, which implies good conditions for seismic facies discrimination in this particular geologic scenario.

The Bayesian facies classification was based on several PDFs as functions of Ip and VP/VS for each facies to address the relationships between the seismic parameters and reservoir properties (Doyen, 2007). These are, however, only based on training data from Margretheholm-1A, which may inadequately represent other possible geologic scenarios in the prospect area, located at a distance of 30 km. However, additional petrophysical logs are available from Karlebo-1A, which can be used to further extend the training data set beyond the well control using statistical rock physics (Mavko and Mukerji, 1998; Avseth et al., 2001; Mukerji et al., 2001; Eidsvik et al., 2004; Doyen, 2007). Figure 5a shows the porosity plotted against the shale volume data from Margretheholm-1A and Karlebo-1A combined, as well as the prior distributions for each lithofacies. Figure 5b shows an augmented data set generated from a Monte Carlo simulation, which randomly samples from the prior distributions. This augmented data set can subsequently serve as the input for the facies-dependent rock-physics models that computes the seismic parameters, such as the Ip versus VP/VS crossplot in Figure 5c, where the PDFs represent well data extrapolated using the statistical rock physics for each lithofacies. We note the occurrence of high-probability areas for the nonreservoir PDF, which represents various shales from the cap rock and the intrareservoir in the Margretheholm-1A well. Finally, these PDFs in the Ip versus VP/VS domain can be used to analyze the seismic inversion results.

For statistical rock physics, different facies require tailored rock-physics models to link the desirable reservoir properties with their seismic parameters. Obtaining the appropriate rock physics model candidates significantly relies on the consolidation degree of a specific rock. Regional basin modeling and core analysis (Japsen et al., 2007; Weibel et al., 2017) have indicated the occurrence of small volumes of quartz cement in the Gassum sandstone formation that formed due to chemical compaction before subsequent uplift events. The Gassum Formation has an estimated burial depth of >2250 m from wells further southwest at Zealand (Weibel et al., 2017), at which temperatures were close to 70°C and the onset of quartz cementation (Avseth et al., 2010). From testing a range of various rock-physics models, we finally selected the pressure-insensitive constant cement model (see Appendix A for details), which appears to predict the seismic velocities of the LJRU and Gassum sandstones quite well for the target depth interval. Figure 6a and 6b shows the petrophysical porosity, fluid saturation, and shale volume logs from Margretheholm-1A, which were used as the model input for the P- and S-velocities and density logs in Figure 6c6e. Model accuracy is generally higher within the cleaner sandstone zones (zones 2, 4, and 6), but the model can also provide sufficient velocity predictions in zones that are more shale dominated. To address model uncertainties associated with natural variabilities at a distance from the well locations, model parameters were expressed with normal distributions (Grana, 2014) defined by a mean value and standard deviation to sample the random values generated by the Monte Carlo simulations.

Furthermore, we calibrated a rock-physics model for the Fjerritslev Formation cap rock using a constant-clay model (Avseth et al., 2005). This model is based on Hertz-Mindlin contact theory (Mindlin, 1949), which is normally used to model high-porosity sands. Despite the model do not address the complex microstructure and anisotropy of shales, it has been demonstrated to be useful for modeling sandy shales (Avseth et al., 2005), and we found it to yield sufficient velocity predictions for the Fjerritslev Formation as well. Therefore, this model was used to extend the nonreservoir PDF (see “Nonreservoir model” in Figure 5c).

Seismic inversion

Prior to the seismic inversion, the key seismic processing steps were anomalous amplitude, linear, and coherent noise attenuation; minimum-phase filtering; despiking and predictive deconvolution; datum correction; zero-phase filtering; full Kirchhoff prestack time and depth migration; stacking velocity analysis; normal moveout correction; residual amplitude analysis and compensation; estimation of angle gathers; wide-angle mute; and time alignment. In general, land-seismic acquisition is hampered by the complex and heterogeneous near-surface zone of the earth due to distortion of signals caused by irregular time delays, scattering, and energy absorption, leading to loss of high-frequency signals. In addition, surface waves (e.g., ground roll) can generate coherent surface-related noise. Despite applying different approaches to compensate for these effects (e.g., static correction of time delays), the signal-to-noise ratio (S/N) is often not satisfactory improved (Schmelzbach et al., 2016a).

In particular, the seismic lines in our data set suffer from varying fold coverage due to terrain obstacles and local areas where seismic data are unavailable. Because the seismic data are inverted in the angle-gather domain, the locations of areas with low or no folding vary spatially and temporally between each angle stack, which introduces possible inversion errors and instability. Therefore, an S/N volume was applied to weight the data contribution that was proportional to the fold magnitude in the inversion, as shown in Figure 7 for line 5. The distance between the Karlebo-1A well and the nearest CMP seismic gather is approximately 140 m, located in a low-fold zone, which indicates that the well tie and calibration procedures are more demanding for the seismic inversion. Therefore, the well was relocated further to the northeast (i.e., to the right in the seismic profiles in Figure 7), where a larger seismic fold is present and more representative of the other inverted lines, despite increasing the distance to approximately 600 m.

In this study, we performed a simultaneous seismic AVO inversion based on 2D seismic offset gathers stacked into angle bands of 3° to derive the absolute acoustic impedance Ip, the P-to-S velocity ratio VP/VS, and the density (refer to Appendix B for more details on this method). The background models were primarily developed based on the Karlebo-1A well logs as the seismic inversion results benefit from a comparison with and verification by the nearest well data. In addition, the constrained low-frequency seismic content rendered the background models even more dependent on the adjacent well-log measurements. Because the S-velocity and density logs are unavailable for the Karlebo-1A well due to technical issues, we empirically derived these data using the Margretheholm-1A well as a calibration point. In the wavelet estimation, we extracted the statistical wavelets from the seismic frequency spectra (Edgar and Baan, 2011) that were convolved with the angle-dependent seismic reflectivities. The phases of the wavelets were estimated by correlating the inversion results to the well-log data.

Figure 8 displays, from left to right, the Karlebo-1A well logs, background models, and the seismic inversion results projected onto the relocated Karlebo-1A well. Due to a low maximum angle of incidence of up to 27°, reliable density inversions are problematic, and the VP/VS inversions may become unstable and uncertain, as the noisy VP/VS results demonstrate in Figure 8. Normally, angles of up to 40° can give good VP/VS inversion results when using a three-term AVO approximation (e.g., Aki and Richards, 1980) but, in this case, the high-velocity contrast at the base of the chalk produces refraction effects beyond an estimated 35° angle, which hampers the primary signals from the reservoir boundaries. Despite these issues and data limitations, the inverted Ip and VP/VS are reasonably balanced and in good agreement with the elastic logs from Karlebo-1A. Immediately below the Base Chalk formation, however, we observed several extreme VP/VS values that are associated with local extreme values in the seismic data input. Nevertheless, we suggest that the acoustic impedance and VP/VS results within the target LJRU and Gassum Formation are reasonably well balanced.

Reservoir characterization

In geothermal exploration and derisking, porosity and permeability are the two key parameters used to characterize a reservoir. These can serve as input to reservoir simulations to evaluate the potential heat extraction, storage, and production lifetime of geothermal reservoirs. Furthermore, the reservoir performance highly depends on the lithologic composition of the reservoir, for example, higher shale volumes within a sandstone unit will typically reduce the permeability and reservoir quality. However, we are unable to infer all the reservoir parameters from the seismic data. For example, the permeability does not have a direct influence on rock stiffness and is, therefore, derived from the porosity predictions in this study (see Appendix C for details). Furthermore, the pore fluid temperature, pressure, and chemical composition are also interesting parameters in a geothermal context. However, based on the rock-physics model calibrated to the LJRU and Gassum Formation, sensitivity studies showed relatively small effects on the seismic response when varying these parameters. Hence, this study is limited to predicting the spatial reservoir quality distributions, in terms of facies, porosity, permeability, and shale volume, to identify the most promising geothermal sites.

We combine the results from the statistical rock physics and seismic inversion to quantitatively characterize the geothermal target reservoirs in a Bayesian probabilistic framework (refer to Appendix C for a more detailed description of this approach). Before we consider the seismic inversion data, we must test the reliability of our reservoir predictions that were conditioned by the acoustic impedance (Ip) and P-to-S velocity ratio (VP/VS) measurements from Margretheholm-1A. Figure 9 shows, from left to right, the Ip and VP/VS logs; the posterior probability of each facies (or facies probability); the most probable (or predicted) and reference facies profiles; and the posterior distributions of porosity, horizontal reservoir permeability, and shale volume represented by the weighted mean (the red curve) and a 95% confidence interval (the light-blue area), with the petrophysical logs (the black dots) plotted at the top. The net-to-gross volumes were used as prior information in the facies classification.

To evaluate the facies classification accuracy in Figure 9d, which was conditioned by the reference facies in Figure 9e, Figure 10 presents a summary of the results in a confusion matrix (Grana et al., 2017). The diagonal and off-diagonal elements represent the classification success and failure scores, respectively. For data samples that represent the nonreservoir, shaly sandstone, and clean sandstone facies, 82.5%, 40.6%, and 66.9% were correctly classified, respectively. The highest misclassifications were associated with the shaly sandstone facies. For example, 41.4% and 18% of the data classified as shaly sandstone were actually from the nonreservoir and clean sandstone facies, respectively. Moreover, 21.5% of the shaly sandstone was misclassified as clean sandstone. These results demonstrate the nonuniqueness of the problem and imply that there is a certain risk of disorientating the various facies.

Furthermore, there was generally good agreement between the posterior distributions of porosity, permeability, and shale volume and the petrophysical logs in Figure 9f9h. Wider confidence intervals imply less constrained and more uncertain solutions. Within the clean and shaly sandstone facies, the weighted mean porosity tended to be slightly underestimated but the petrophysical logs still plotted within the confidence intervals. We observed an identical trend for the permeability predictions because these were directly derived from the porosity solutions. Higher shale-volume predictions correlate well with lower permeability because an abundance of clay will typically choke the pore throats and reduce fluid mobility. The predictions within the nonreservoir facies also match reasonably well with the petrophysical logs, despite the application of a less rigorous rock-physics model for the shales.

Figure 11 shows the sections of the (Figure 11a) full-stack data, (Figure 11b) inverted Ip, and (Figure 11c) inverted VP/VS along line 5 in a southwest (SW)–northwest direction (see Figure 1). The white stippled, gray curve is the projection of the Karlebo-1A well onto the sections. Based on these seismic inversion results, Figure 12 shows, point by point, the categorically estimated (Figure 12a) facies classification with a gray seismic backdrop and the posterior probabilities of (Figure 12b) clean sandstone, (Figure 12c) shaly sandstone, and (Figure 12d) nonreservoir facies, which cover the target and cap-rock intervals. The results reveal several high-probability lateral distributions of clean and shaly sandstones that are mainly within the LJRU and Gassum Formation. Due to the similar geologic characteristics and elastic properties of clean and shaly sandstones, the nonreservoir facies profile provides a possibly improved sandstone indicator represented by low posterior probabilities. Based on our predictions, the Karlebo-1A well is in a good position to encounter the high-probability sandstone units along this seismic profile. The clean and shaly sandstone facies were modeled using pore fluid properties of in situ brine, but we also investigated other possible pore fluids, such as oil and gas, to better capture the uncertainties. However, the probabilities for brine pore fluids were higher within our target geothermal reservoir, which is also reasonable because there were not hydrocarbons discovered in the Karlebo-1A well.

Furthermore, Figure 13 shows the continuous reservoir parameters using the weighted means of porosity, horizontal reservoir permeability, and shale volume in the first column and, correspondingly, for the weighted standard deviation in the second column. The facies classification results from Figure 12a were used to select, point-by-point, an associated rock-physics model. By considering several of the more promising clean sandstone units, the estimations typically ranged between 20% and 25% for porosity, 400 and 700 mD for permeability, and 16% and 20% for shale volume. These reservoir predictions should support a sufficient geothermal production performance. Furthermore, the corresponding weighted standard deviation of the posterior distributions has a maximum of approximately 3% porosity, 250 mD permeability, and 15% shale volume. Based on this, we can obtain similar uncertainty trends for the 95% confidence intervals in the well-log data (see Figure 9). More uncertainty is associated with the permeability compared with the shale volume, which, in turn, is more uncertain than the porosity. Nevertheless, we predict that the most promising sandstone beds are located within the LJRU, as well as in the upper part of the Gassum Formation. Although the clean sandstone probability slightly decreased when moving away from the Karlebo-1A location, the reservoir properties show little lateral variability within the predicted sandstone bodies (see Figure 12b).

To evaluate the prediction accuracy, we compared the results from the seismic inversion, facies classification, and estimated reservoir properties along the Karlebo-1A projection with the upscaled logs. Figure 14a and 14b shows the seismically inverted Ip and VP/VS values extracted from Figure 11b and 11c, where the VP/VS ratio is characterized, as expected, by several more unstable results than Ip. Figure 14f and 14g shows the smoothed porosity and shale volume solutions extracted from Figure 13a and 13c. These predictions correlate strongly with the seismic inversion; therefore, we observe corresponding misalignments with the logs. The seismic inversion and reservoir parameter predictions, however, generally match the logs.

Discussion

This study demonstrates how seismic reservoir characterization can provide supplementary quantitative facies and reservoir information in a geothermal exploration setting. However, exploiting the full potential of this approach relies significantly on data quality and availability, which are quite restricted in our study area. For example, the missing density and S-velocity logs for the Karlebo-1A well yields a reduced number of locally independent data parameters, which reinforce the nonuniqueness of the inverse problem. Moreover, the quality of the seismic data suffers due to low-fold zones, constrained low-frequency content, remnant multiples, and low maximum offsets. Therefore, our results have a significant uplift potential when we have access to a complete set of well logs, which can provide and restore the low frequencies and high offset data during seismic acquisition and processing. Based on this, we demonstrate that it is possible to derive reasonable reservoir predictions, despite a restricted data foundation.

To perform assessments of our reservoir characterization, the confusion matrix in Figure 10 indicates that there is a nearly 30% probability for an incorrect facies prediction of the data samples in Margretheholm-1A. This is due to the inherent overlap of seismic properties (i.e., the acoustic impedance and P-to-S velocity ratio) for the various facies, which increases the risk of confusing one facies with another. In most cases, we are able to distinguish between nonreservoir and reservoir sandstones, which is encouraging, but this also implies that there are particular challenges associated with classifying the shaly sandstone. In the Margretheholm-1A well, the reservoir unit partly consists of intrashale sequences with elastic properties that differ from the cap-rock interval but have more similar properties to the shaly sandstone (see Figure 5c). These nonunique characteristics yielded the 41.4% misclassification of the actual nonreservoir as shaly sandstone. Therefore, combining two different types of shales into one single nonreservoir PDF is, perhaps, an oversimplification of the facies classification. However, at a seismic resolution, a less detailed classification scheme allows us to outline the target reservoir sandstone units (Mukerji et al., 2001; AlKawai et al., 2018). Alternatively, we can further simplify the facies classification scheme by combining the clean and shaly sandstone PDFs into a single-reservoir PDF. On the other hand, this would prohibit us from distinguishing the variations in the reservoir quality.

A second assessment of our results is based on the Karlebo-1A data, as shown in Figure 14, where we observe a varying mismatch between our estimates and the well-log measurements. There are several possible explanations for this, but perhaps the most relevant explanation is that the position of the Karlebo-1A well is at a distance of approximately 600 m to the CMP seismic gather that was used for the well tie. This CMP seismic gather was selected because the nearest gather is located in a low-fold zone (see Figure 7a), which provided an insufficient well tie for the seismic inversion. Consequently, any direct comparison between a seismic trace and the well logs is uncertain, especially as the subsurface is characterized as densely faulted and laterally variable. In addition, an imperfect time-to-depth relationship and a deviated well path through the target reservoirs may contribute to this mismatch. For example, we predicted that several sandstone units are below the Base LJRU, which appears to not be present in the Karlebo-1A logs. Furthermore, the underpredicted shale volumes within the more shaly sections are, presumably and in part, due to the slightly lower shale volumes observed in the Margretheholm-1A well, which were used as training data. This is a reminder of the possible consequences associated with the inclusion of training data from wells at a distance, which may be characterized by slight deviations in the geologic scenario. However, implementing spatial constraints is generally problematic and even more so with a low local well coverage. Also, the lithology predictions are highly dependent on the inverted P-to-S velocity ratio accuracy, which can be, as previously mentioned, unstable due to limited offset and absent S-velocity measurements in Karlebo-1A.

The facies predictions show small lateral variabilities when moving away from the Karlebo-1A position. Figure 12b shows a decreasing probability for the clean sandstone toward the southwest along line 5 in the LJRU, which may imply that the reservoir quality decays away from the Karlebo-1A area. At the same time, we observed that these predictions were distributed in a more discontinuous and patchy manner, which is possibly related to seismic interference or resolution issues. Also, the nonparametric, Kernel-density-estimated PDFs, used to link reservoir and seismic properties (see Figure 5c), were generated using a calibration procedure from the well-log data (Grana, 2018). Obtaining optimal parameters to control the PDFs can be challenging for non-Gaussian distributed data, such that we finally selected a configuration that yielded predictions that were consistent with the expected geology. Meanwhile, during PDF calibration, we observed that different configurations influenced the predictions in terms of varying spatial continuity, resolution, and noise. With such instabilities, there was a significant reduction in our prediction confidence at increasing distance from the Karlebo-1A location. Although various statistical estimators can be selected to provide a reasonable representation of the reservoir predictions, we applied a weighted mean and standard deviation (see Figure 13). Previous studies have demonstrated, for instance, the use of the maximum a posteriori (MAP) estimator (Bachrach, 2006) and various percentiles of the posterior PDFs (Grana, 2018). Nevertheless, we obtained similar results when testing various estimators.

Conclusion

In this study, we demonstrated a field example of seismic reservoir characterization in a low-enthalpy geothermal exploration setting from an area to the north of Copenhagen, onshore Denmark. Facies and reservoir property predictions were derived by integrating seismic data, well logs, and available regional data. The results demonstrate that several porous and clean water-bearing sandstones are potential high-quality geothermal reservoirs within two target layers, that is, the LJRU and Gassum Formation. In general, the LJRU exhibits a more promising reservoir quality in terms of high porosity, permeability, and low shale content compared with the Gassum Formation. However, with the limited availability and nonoptimal quality of the local data, we were unable to derive a rigorous and firm reservoir prognosis and our final results are primarily indicative of possible target layers. With an improved data foundation, geothermal reservoir seismic characterization is demonstrably useful to optimize site locations, which can be the tipping point between success and failure in geothermal exploitation.

Acknowledgments

This contribution is a part of the project GEOTHERM (Geothermal energy from sedimentary reservoirs — Removing obstacles for large-scale utilization) funded by the Innovation Fund Denmark (IFD) (grant no. 6154-00011B).

Data and materials availability

Data associated with this research are available and can be obtained by contacting the corresponding author.

Rock-physics modeling

For isotropic rocks, the effective bulk modulus K and shear modulus μ are related to the seismically inverted acoustic impedance Ip and the ratio of the P-wave velocity VP to S-wave velocity VS (i.e., VP/VS) based on the following simple expressions:  
VP=K+43μρ,
(A-1)
 
VS=μρ,
(A-2)
 
Ip=VPρ.
(A-3)
The effective density ρ is estimated with the following equation:  
ρ=(1ϕ)ρm+ϕρf,
(A-4)
where ϕ is the porosity and ρm and ρf are the mineral and fluid densities, respectively.
The elastic modulus for a dry rock in a high-porosity regime is estimated based on the contact cement theory (CCT) (Dvorkin and Nur, 1996), which is valid for cemented sediments:  
KCCT=16n(1ϕ0)McSn,
(A-5)
 
μCCT=35KCCT+320n(1ϕ0)μcSτ,
(A-6)
 
Mc=ρcVpc2,
(A-7)
 
μc=ρcVsc2,
(A-8)
where ϕ0 and n are the critical porosity and coordination number, respectively; ρc, VPc, VSc, and Mc are the density, P-velocity, S-velocity, and compressional modulus of the cement component, respectively; and Sn and Sτ are the normal and shear stiffness coefficients, respectively. We model a constant 3% porosity reduction from ϕ0 to ϕCCT due to the quartz cement in the CCT. We then interpolate, for a given porosity ϕ, between the high-porosity end member (i.e., ϕCCT=0.37) and the zero-porosity mineral point using a lower Hashin-Shtrikman-Walpole (HSW–) (Walpole, 1966a, 1966b) to model the bulk KHSW and shear μHSW moduli (Mavko et al., 2009) with the following equations:  
KHSW=(ϕϕCCTKCCT+43μCCT+1ϕϕCCTK+43μCCT)143μCCT,
(A-9)
 
μHSW=(ϕϕCCTμCCT+z+1ϕϕCCTμ+z)1z,
(A-10)
where K and μ are the bulk and shear modulus of the solid mineral, respectively, and z is defined as follows:  
z=μCCT6(9KCCT+8μCCTKCCT+2μCCT).
(A-11)
The pore fluid effect is modeled using the equations reported in Gassmann (1951), whereas the equations reported in Voigt (1928) approximate a patchy fluid mixture (Mavko et al., 2009). Table A-1 lists the model parameters defined by the mean values M and standard deviations σ. The calibrated high fluid densities are possibly due to uncorrected mud invasion that affected the density measurements. Furthermore, high brine salinity consequently yields a correspondingly high bulk modulus. Table A-2 summarizes the remaining model parameters.

Simultaneous AVO seismic inversion

A simultaneous AVO seismic inversion converts the prestack seismic data into absolute elastic parameters, such as acoustic impedance, seismic velocities, and density (Mallick and Adhikari, 2015). In this study, the AVO inversion uses the three-term linearized approximation of the Zoeppritz AVO equations reported in Aki and Richards (1980) to model the P-wave reflection coefficient Rpp as a function of the angle of incidence θ:  
Rpp(θ)12(ΔVPVP+Δρρ)[12ΔVPVP2VS2VP2(2ΔVSVS+Δρρ)]sin2θ+12ΔVPVP(tan2θsin2θ),
(B-1)
where θ=1/2(θ1+θ2)θ1, ΔVP=VP2VP1, VP=1/2(VP1+VP2), ΔVS=VS2VS1, VS=1/2(VS1+VS2), Δρ=ρ2ρ1, and ρ=1/2(ρ1+ρ2); and VP, VS, and ρ, respectively, denotes P- and S-velocities and density, and indices 1 and 2 denote the top and base layer, respectively. Equation B-1 yields robust AVO approximations to the Zoeppritz equations up to approximately 40°. By sequentially updating a low-frequency background model of velocities and densities obtained from the upscaled well-log data and seismic interval velocities, we obtain an adequate match between the forward modeled and real seismic data. Hence, we can systematically approach velocity and density models of the subsurface with higher confidence.

Facies and reservoir property estimation

In this study, we estimated the facies and reservoir properties in a Bayesian framework. The reservoir properties were estimated from probability distributions, which were explicitly expressed using PDFs to conform to the nonlinear relationships between the reservoir and seismic properties (Grana, 2018). This problem can be formulated as a Bayesian inference problem:  
σ(m)=cp(m)L(dobsg(m)),
(C-1)
where m denotes the subsurface reservoir properties and facies and c is a normalization constant. Data for m prior to including seismic observations are expressed by the prior probability distribution p, which is normally obtained from adjacent well data. The model likelihood function L was then used to implement seismic data, which measure, in terms of probability, the mismatch between the modeled g(m) and seismically inverted data, dobs. In this case, g(m) and dobs refer to the PDFs as a function of the augmented Ip and VP/VS training set and seismic inversion data, respectively. We use a nonparametric kernel density estimation technique to calculate the PDFs (Silverman, 1986; Grana, 2018).
Based on the posterior distributions, the most probable facies was estimated from the MAP probability (Doyen, 2007). Furthermore, we calculated the weighted mean and standard deviation from the most likely facies posterior distribution to quantify the expected values and uncertainty ranges for the porosity and shale volume solutions. The weighted mean for the porosity ϕF in facies F is expressed by the following equation:  
ϕ¯F=k=1LP(ϕFk|Ip,VPVS)ϕFkk=1LP(ϕFk|Ip,VPVS),
(C-2)
where k is a modeled porosity sample and P(ϕk|Ip,VP/VS) is the posterior probability for a set of Ip and VP/VS data. The corresponding weighted standard deviation is expressed by  
σ¯ϕF=k=1LP(ϕFk|Ip,VPVS)(ϕFkϕ¯F)2k=1LP(ϕFkIp,VPVS).
(C-3)

Permeability estimations

To estimate the horizontal reservoir permeability κres, we use an empirical permeability-porosity model κL as proposed by Kristensen et al. (2016):  
κL=c×(4.43×104×ϕ)4.36,
(C-4)
where ϕ is the predicted porosity and c is a constant controlled by variations in the local permeability calibrated to 0.85 in this study. To convert the relationship in equation C-4 from laboratory measurements of horizontally orientated core plugs to field-scale reservoir permeability κres, κL is multiplied by an upscaling factor of 1.25. For the corresponding vertical permeability, core analyses suggest that a vertical-to-horizontal permeability ratio of 0.3 is appropriate.

graphic
Kenneth Bredesen received a Ph.D. (2017) in geophysics from the University of Bergen. He is a geophysicist at the Geological Survey of Denmark and Greenland (GEUS). He has previously worked as a QI geophysicist at Spike Exploration and as a postdoctoral researcher at Aarhus University. His core competencies are quantitative seismic interpretation, rock-physics modeling, seismic attribute, and amplitude analysis.

graphic
Esben Dalgaard received a Ph.D. (2014) in geophysics from Aarhus University. During his Ph.D., he worked on developing and implementing analogue and digital processing techniques for geophysical instruments. He is a senior geophysicist at Qeye in Copenhagen. At Qeye, he primarily works with seismic inversion and rock-physics solutions for quantitative interpretation on data from around the world.

graphic
Anders Mathiesen is a senior advisor at GEUS. He has previously worked with geophysical geology including research, advisory, and consult work within hydrocarbon, CO2 storage, and geothermal utilization. His core competencies are geophysical and geological data integration, analysis, implementation and mapping, modeling, and assessment based on various data (seismic, cores, logs, etc.).

graphic
Rasmus Rasmussen is a senior advisor at GEUS. He has been working mainly with seismic data processing offshore and onshore and also with seismic acquisition, seismic inversion, and AVO processing. In addition to this, he has been working with seismic forward modeling for geologic model verification. The work has been related to hydrocarbon exploration projects, geothermal projects, and CO2 storage projects.

graphic
Niels Balling received a Ph.D. (1982) in geophysics from Aarhus University and was offered the doctor of science degree in 2013. He is an associate professor of geophysics in the Department of Geoscience, Aarhus University. He is currently sector leader of ‘Deep Earth Systems’. He has core competences within several geoscience disciplines including controlled and natural source seismology and thermal modeling, all in relation to basic crustal and lithospheric studies as well as geothermal energy resources.

Freely available online through the SEG open-access option.