Deep learning for end-to-end subsurface modeling and interpretation: An example from the Groningen gas field

While playing a key role in the cycle of reservoir field exploration and development, reliable interpretation of subsurface geology is a complicated, time-consuming process that usually requires interdisciplinary collaboration between geophysicists, geologists, petrophysicists, and petroleum engineers. The workflow traditionally follows a manual and labor-intensive approach that aims to integrate various subsurface data, including seismic, well log, and production, through a set of parallel or consecutive subtasks, such as data conditioning and processing, structure mapping, stratigraphy delineation, property modeling, etc. For automating these various tasks, in addition to the large number of traditional workflows and algorithms, more tasks have been developed based on the recent success of deep learning, particularly the convolutional neural networks (CNNs) used in language, image, and video processing. These applications include log interpretation (e.g., Tschannen et al., 2017; Wu et al., 2018; Zhu et al., 2018; Kaul et al., 2021a), fault detection (e.g., Huang et al., 2017; Di et al., 2018a; Guitton, 2018; Xiong et al., 2018; Zhao and Mukhopadhyay, 2018; Wu et al., 2019a; Zhao, 2019), salt body isolation (e.g., Di et al., 2018b; Shi et al., 2018; Waldeland et al., 2018; Kaul et al., 2021b), horizon interpretation and stratigraphy mapping (e.g., Li et al., 2019; Peters et al., 2019; Wu et al., 2019b; Di et al., 2020a, 2021; Li and Abubakar, 2020; Zhao et al., 2021a), property estimation (e.g., Alfarraj and AlRegib, 2019; Das et al., 2019; Mustafa et al., 2019, 2020; Wang et al., 2019; Cai et al., 2020; Di et al., 2022), and so on. To the best of our knowledge, each application is often performed individually, and little effort has been made toward an integrated workflow of all components essential to subsurface interpretation. For example, Wang et al. (2018), while summarizing these relevant interpretation applications, do not attempt to integrate them into an end-to-end workflow. Di et al. (2020b) demonstrate how the CNN greatly assists fault and stratigraphy interpretation in the Opunake field, but they do not cover other aspects for comprehensive subsurface interpretation.

In this paper, we present such a deep-learning-powered end-to-end workflow capable of automatically conditioning well logs, detecting faults and salts, estimating relative geologic time (RGT), and generating and validating rock property models over the Groningen gas field.

The Groningen gas field, discovered in 1959 and in operation since 1963, is located to the northeast of the Netherlands in the southern North Sea. It produces natural gas from the Permian Rotliegend Group of an approximately 50–225 m thick Slochteren sandstone, which is covered by thick Zechstein salt layers, distorted by a complex fault system, and composed of mainly Aeolian and fluvial sandstone interbedded with silt and shale of excellent reservoir properties (porosity approximately 15%–20%; permeability approximately 0.1–3000 mD) (Figure 1a) (Evans et al., 2003; van Thienen-Visser and Breunese, 2015; van Thienen-Visser et al., 2015). The released data set (courtesy of Utrecht University) consists of a 3D seismic poststack amplitude volume; more than 400 wells of logs including sonic, density, gamma ray (GR), porosity logs, and more; fault and horizon interpretation of multiple sub- and supra-Zechstein surfaces; and a reservoir model of the Rotliegend reservoir formation (Figures 1b and 1c). The seismic covers an area of approximately 2700 km² with 1961 inlines and 2181 crosslines and is well imaged particularly above the Zechstein salt layer. All wells have caliper and GR logs. The density log is available in most of the wells but exists in short depth ranges. The sonic log is relatively more complete than density in each well but only exists in 160 wells, and the included porosity log is limited to only the Rotliegend reservoir formation where the log is present.

Starting with the available seismic and well logs, the proposed end-to-end workflow for deep-learning-assisted subsurface interpretation of the Groningen field consists of the following six components (Figure 2):

Log reconstruction. The first workflow component aims at the reconstruction of data from the missing well logs to make them more complete. Figure 3 illustrates the workflow for CNN-based reconstruction of density and sonic logs in the Groningen data set, which uses an encoder-decoder architecture for capturing the relationship between the GR, density, and sonic logs (Chen et al., 2019). For CNN training and validation, we first split the 153 wells with all three logs into 130 training wells and 23 validation wells.

The workflow is first performed for reconstructing the missing density log from GR and sonic logs. Figure 4a displays the results on four of the 23 validation wells, including AMR-1, ROD-201A, SWO-1, and TBR-4. Note that AMR-1 and ROD-201A are in proximity of other wells, while SWO-1 and TBR-4 are located where there is no other well in the neighborhood. A good agreement is observed between measured and predicted logs for most depths. For regions such as near 3000 m in well SWO-1 and near 3150 m in well ROD-201A, the mismatch between predicted and logged density agrees with alterations in the caliper log in these wells and is likely caused by a washout effect; the predicted density is closer to formation density in this interval. Similarly, the same workflow is applied for sonic log reconstruction from GR and density logs, with the results at the same four validation wells shown in Figure 4b. The final error in the test wells is comparable with the tool's uncertainty in most of the depths, and that quality is maintained even in the sand reservoir section, presenting a higher vertical heterogeneity compared with the seal, for example in well TBR-4. For better evaluating the predicted log, the depth-by-depth uncertainty helps identify whether the input logs provide enough information to predict the remaining log and how confident the model is in predicting the missing information.

After validation and improvement, the trained CNN model is applied to reconstruct the missing logs in the entire data set, resulting in a complete set of GR, density, and sonic logs in all of the wells.

Salt body identification. The Groningen field consists of several fault structures that are a result of salt tectonics in the area, so it is important to understand the placement of salt pillows, also known as Zechstein salt (Logeman, 2017). For identifying the salt bodies throughout this 3D seismic, our workflow is CNN-based image segmentation, which is in the architecture of U-net with dilated convolution (Figure 5) (Kaul et al., 2021b). For training the network, we identify and annotate the top and base of the Zechstein salt in 13 inline and 15 crossline sections, which amount to 0.67% of the total seismic available. More specifically, the base of salt is of a strong impedance contrast due to the direction interaction between the salt and the reservoir sandstone below, whereas the top of salt has weaker impedance contrast due to the mixture of salt and sediments, as shown in one of the training sections (inline #8100) (Figure 8). After training the CNN using the 28 annotated sections, we run the trained CNN on every inline and crossline section, which provides two salt body probability volumes. Both are further consolidated as the root-mean-squared value at each sample to produce a smoother salt body prediction and to minimize false positive predictions that are present in one direction but not in the other.

Figure 6 shows both top and bottom views of the 3D salt body predicted by the CNN workflow, which clearly highlights the presence of multiple salt pillows and the distortion by the complex fault network below the Zechstein salt.

Fault detection. Given the strong connection between the Zechstein salt, the Rotliegend reservoir, and the overburden faults, it is important to interpret the supra-Zechstein fault system. Figure 7 illustrates our CNN-assisted fault detection workflow, which consists of two subnetworks and is unique in the industry (Li and Abubakar, 2019). Specifically, for a given 3D seismic volume, interpreters manually interpret fault labels on several inlines and crosslines, which are used to train the first network (CNN#1). Then the trained CNN#1 is applied on every inline and crossline to produce two fault predictions. These predictions are then fused together and fed into the secondary network (CNN#2) to generate the final prediction volume. The CNN#2 has been trained to process time slices to enhance feature connectivity and reduce overall noise. This dual-network solution has been tested capable of producing fault detection of higher quality based on a very small amount of manual interpretations (labels). More technical details about the CNN implementation for seismic fault detection can be found in Li and Abubakar (2019).

For detecting the supra-Zechstein fault network in the Groningen field, we train CNN#1 using only five inlines and four crosslines (Figure 8a), which efficiently images the complex fault system throughout the entire Groningen field (Figure 8b).

Relative geologic time estimation. An RGT volume is a seismic attribute volume in which each sample is assigned a relative time (or geologic age) corresponding to the seismic amplitude at the same sample. Such a volume successfully preserves the spatial relations and temporal ordering between all structures and stratigraphic units in the survey and is proven to be a strong constraint to CNN-based seismic interpretation and property inversion. Therefore, RGT estimation is an essential component of the proposed end-to-end subsurface interpretation. Figure 9 illustrates the CNN-assisted workflow for self-supervised RGT estimation, which uses the Seismic FlowNet (Li and Abubakar, 2020) architecture to estimate the relative displacement between adjacent seismic sections. Moreover, to avoid dependence on manual annotation, the training data are randomly extracted from the target 3D seismic and deliberately warped by predefined 2D flow fields. In this way, no user input label is required. The corresponding RGT volume is obtained by extracting a dense set of horizons using the estimated flow fields. More technical details about the CNN implementation for Seismic FlowNet can be found in Li and Abubakar (2020).

Figure 10 shows the densely picked horizons and corresponding RGT volume of the Groningen field, respectively, both of which outline the structures in this field well, including the Zechstein salt, the overlying deformation, and the underlying Rotliegend reservoir.

Property estimation. With the provided seismic and well logs as well as the generated RGT and structure volumes, the next component in our subsurface interpretation workflow (Figure 2) is property estimation or inversion. The corresponding workflow is shown in Figure 11. It consists of two steps: unsupervised learning for seismic feature engineering and supervised learning for seismic-well integration. Each step is solved using a CNN. The first step utilizes an unsupervised learning network for extracting the regional features present in the target area by learning the input seismic (amplitude and/or attribute), fault prediction, and RGT volumes, while the second CNN builds the optimal nonlinear mapping between the four target properties and the seismic (amplitude) patterns at the well locations. Both components are connected by embedding the first network into the second one. More technical details about the semisupervised property inversion can be found in Di et al. (2022). For validating the machine predictions, 22 of the available wells are reserved from performance evaluation, while the rest are used for training.

Figure 12a compares the machine predictions and the actual log measurements at four of the 22 validation wells, including AMR-1, LRM-7, SWO-1, and TBR-4. The orange band represents the uncertainty associated with the machine prediction, which is estimated using the dropout-based scheme (Zhao and Chen, 2020). Furthermore, we show the predicted sonic and density over inline #8813 of the validation well SWO-1 in Figures 12b and 12c, respectively, both of which are in good match with the seismic patterns and of high lateral continuity throughout the survey. From these results, not only is a good match observed between all four predicted properties and the actual logs, but also the machine successfully propagates its learning throughout the entire section and provides consistent predictions based on the seismic patterns, including the Zechstein salt and moreover the Rotliegend reservoir layer.

Property estimation result verification. To evaluate the property estimation results in the previous component, we next use a machine learning (ML) model trained with reversed data-label pairs, which take rock properties as input and seismic data as target, to qualitatively evaluate the property estimation result (Zhao et al., 2021b). Such evaluation works under the assumption that after training at well locations using the ground truth impedance-seismic pairs and applying to the impedance volume predicted from the previous ML model, the quality of the simulated seismic data qualitatively indicates the quality of the ML-based impedance prediction. Figure 13 summarizes the relation between a property estimation model and the corresponding backward model. On the Groningen data set, we use acoustic impedance as the property for the backward ML model, which is derived from the density and sonic logs. The backward ML model shares largely the same architecture as used in the property estimation tasks, with the exception that it now operates on single traces only, given that the input to the ML model is only available along wellbores. The average trace-wise correlation between the original and predicted seismic amplitude at validation well locations is about 0.62.

Figures 14a and 14b show the original seismic amplitude and the corresponding acoustic impedance converted from the predicted density and sonic along inline #8813, while Figure 14c shows the simulated seismic data from the backward ML model. Compared with the original seismic input, Figure 14a, the backward ML model provides only moderate reconstruction to the original data, possibly due to the error within the input impedance volume propagating to the predicted seismic. Moreover, the amount of training data is very small (only at well locations) and has incomplete coverage vertically. To quantify the misfit between original seismic and the predicted seismic, we compute the correlation coefficient within a 51-sample window centered at each sample. In Figure 14d, we show the predicted impedance volume modulated by the computed correlation, in which the brighter impedance corresponds to samples of higher correlation between the original and backward predicted seismic, and vice versa. We observe that the postsalt regions with less complex structures in general have higher correlation, whereas regions with high deformation, within salt, and beyond well penetration consist of much lower correlation. Qualitatively, we can use such correlation to infer the quality of the predicted impedance cube where the ground truth impedance is not accessible. Such backward ML prediction can be repeated if an updated impedance cube becomes available.

An end-to-end workflow has been presented for automated subsurface modeling and interpretation from seismic and well logs. It successfully integrates multiple key interpretation tasks such as data conditioning and processing, structure interpretation, and property estimation and verification through deep learning algorithms. While successfully tested on the Groningen gas field, it could be readily applied to other similar fields, and it could be extended with more necessary components, such as log outlier detection and seismic conditioning, as required by a specific study area.

The authors would like to thank Nederlandse Aardolie Maatschappij (NAM) and Utrecht University for providing the Groningen data set to the public under the Creative Commons Attribution 4.0 license.

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