Wavelet transform‐based seismic facies classification and modelling: application to a geothermal target horizon in the NE German Basin

At the geothermal test site near Groß Schönebeck (NE German Basin), a new 3D seismic reflection survey was conducted to study geothermal target layers at around 4 km depth and 150°C. We present a workflow for seismic facies classification and modelling which is applied to a prospective sandstone horizon within the Rotliegend formation. Signal attributes are calculated along the horizon using the continuous Morlet wavelet transform. We use a short mother wavelet to allow for the temporal resolution of the relatively short reflection signals to be analysed. Time‐frequency domain data patterns form the input of a neural network clustering using self‐organizing maps. Neural model patterns are adopted during iterative learning to simulate the information inherent in the input data. After training we determine a gradient function across the self‐organizing maps and apply an image processing technique called watershed segmentation. The result is a pattern clustering based on similarities in wavelet transform characteristics. Three different types of wavelet transform patterns were found for the sandstone horizon. We apply seismic waveform modelling to improve the understanding of the classification results. The modelling tests indicate that thickness variations have a much stronger influence on the wavelet transform response of the sandstone horizon compared with reasonable variations of seismic attenuation. In our interpretation, the assumed thickness variations could be a result of variable paleo‐topography during deposition of predominantly fluvial sediments. A distinct seismic facies distribution is interpreted as a system of thicker paleo‐channels deposited within a deepened landscape. The results provide constraints for the ongoing development of the geothermal test site.

. (b) Location of vibroseis sources (red stars), explosion sources (black stars), receivers (grey circles) and resulting nominal fold distribution. results from a seismic facies classification and modelling applied to this new 3D seismic data. A sandstone horizon within the Rotliegend formation is used as an example to show the potential of our seismic facies interpretation approach for the characterization of geothermal reservoirs.
Seismic facies analysis is an efficient tool to study variations of signal characteristics within sub-regions of seismic reflection data. Generally, variations in seismic facies are assumed to reflect geological facies variations to be mapped for reservoir characterization studies. The seismic facies is typically determined based on a series of signal attributes measuring properties such as amplitudes, frequency content and their local distribution over time and space (e.g. Taner, Koehler and Sheriff 1979;Marfurt et al. 1998;Chopra and Marfurt 2005;Barnes 2007). A large number of partly redundant signal attributes are available, which even can make choosing the useful parameters difficult in some cases (Barnes 2007). As an alternative to these more classical attributes, the Morlet wavelet transform delivers spectral amplitudes for the full range of measured frequencies (Grossman and Morlet 1984). We consider the wavelet transform attributes as a kind of full waveform representative which can be used for a holistic characterization of the seismic reflection signals. The waveform-derived signal attributes are combined in different ways of multi-attribute analysis. Some approaches try to predict reservoir properties based on correlations between seismic waveform signatures and physical properties determined in boreholes (e.g. Russell et al. 1997;Trappe and Hellmich 2000;Hampson, Schuelke and Quirein 2001). The most common way to combine multi-attribute information is clustering analysis using statistical methods or neural network-based techniques (e.g. Coléou, Poupon and Azbel 2003;de Matos, Osorio and Johann 2007;Pussak et al. 2014;Molino-Minero-Re et al. 2018). Facies types identified by the clustering analysis can be mapped to show their distribution across the geological targets. The interpretation is typically supported by reconciliation with available information from wells (e.g. Aleardi et al. 2015). In addition, the signal characteristics defining the seismic facies types may provide insight into their inherent geological nature (e.g. Pussak et al. 2014). Ultimately, results from classification and mapping of seismic facies can be used to define geometric and petrophysical constraints for the construction of reservoir models. Our workflow includes wavelet transform analysis of reflection signals, classification of wavelet transform-based patterns using self-organizing maps and interpretation of classi-fied patterns supported by seismic modelling. After a short introduction of the geological setting, experiment and data, we describe our method for the generation of data patterns. We also show details of our method of self-organizing maps, which include some modifications to facilitate the definition and visualization of clusters after the learning procedure. For the interpretation of the classification results, we apply seismic modelling based on borehole information. We chose this approach to improve the understanding of the found seismic facies characteristics as a contribution to overcome the kind of black box nature of neural network-based classification.

S T U D Y A R E A A N D D A T A
The study area is located in the NE German Basin (Fig. 1), which represents a sub-basin within the Southern Permian Basin (SPB). The entire SPB system extends from Middle England to North Germany, Poland and the Baltic States. The basin developed in response to thermal relaxation of Early Permian lithospheric thinning and related magmatic destabilization, crustal extension, tectonic subsidence and sedimentation (van Wees et al. 2000). After a phase of Early  Figure 3 (a) Seismic trace from the 3D seismic data. (b) Time-frequency spectrogram from continuous wavelet transform using a short mother wavelet with l = 5/(2π ). Exemplary wavelets are displayed for frequencies of 30, 55 and 80 Hz. (c) Time-frequency spectrum for a longer mother wavelet with l = 10/(2π ). Higher resolution in frequency is related with lower resolution in time and vice versa.
Permian volcanism, basin subsidence and sedimentation started in the Late Permian and lasted until Early Jurassic times. Generally, the North German Basin is considered, besides the northern foreland of the Alps and the Rhine graben, as a prospective region in Germany to make use of geothermal energy. GFZ Potsdam is developing the geothermal test site in Groß Schönebeck as a research platform since 2000. The location had been selected as an exemplary site for geothermal exploitation in the North German Basin based on criteria such as favourable temperature field, adequate porosity and permeability of potential target horizons, vicinity to potential users (Berlin) and access to an existing borehole (Huenges and Hurter 2002). In the E GrSk 3/90 well at the geothermal test site (Fig. 1), Carboniferous rocks are overlain by Rotliegend volcanics, conglomerates, sandstones, siltstones and mudstones (Fig. 2). The porous and permeable sandstones and volcanic rocks of the Rotliegend are considered as potential target horizons for geothermal exploitation based on surrounding temperatures of 150°C. Here we focus on the sandstone formation, as the identification of the volcanic rocks in the seismic data appears less clear. The Rotliegend sandstones formed under semi-arid to arid, sub-equatorial conditions. Topography and climatic conditions were major controlling factors of deposition (Gast et al. 2010). Our study area ( Fig. 1) is located within a transition zone characterized by eolian, fluvial and playa facies (Gast et al. 2010;Fryberger et al. 2011). The major depocentre was located more towards the west (Saline lake in Fig. 1a) while sediment input was presumably delivered from the east where a minor swell separates the NE German Basin from the Polish Trough. The targeted Rotliegend formation is overlain by carbonates, anhydrite and halite of the Zechstein saliniferous successions (Fig. 2) widespread abundant in the SPB. A first geological model was developed by use of borehole information and reinterpretation of sparse two-dimensional (2D) seismic reflection data from former industrial exploration within the region (Moeck, Schandelmeier and Holl 2009). A combined seismic tomography and magnetotelluric experiment within the European Union-funded I-GET project delivered a combined 2D model of seismic velocity and electrical resistivity (Bauer et al. 2010;Muñoz et al. 2010;Bauer, Muñoz and Moeck 2012). A new well (Gt GrSk 4/05) was drilled and used in combination with the E GrSk 3/90 well  for a series of hydraulic stimulation experiments to investigate the feasibility of producing geothermal energy (e.g. Legarth, Huenges and Zimmermann 2005;Blöcher et al. 2016). These studies were accompanied by downhole measurements and core analysis for a petrophysical characteriza-tion of the Rotliegend reservoir (e.g. Trautwein and Huenges 2005).
A new 3D seismic reflection experiment was conducted in 2017 (Stiller et al. 2018;Krawczyk et al. 2019) in order to support further steps in the development of the geothermal research platform in Groß Schönebeck. The experimental layout was centred around the E GrSk 3/90 and Gt GrSk 4/05 wells and covered an area of about 8 by 8 km 2 (Fig. 1b). Source and receiver points with spacings of 50 m were deployed along approximately EW running source lines and NS running receiver lines. Spacing between source and receiver lines was 700 m and 400 m, respectively. A total of 1830 vibroseis points (4 vibrators, linear sweeps of 12-96 Hz) were complemented by 15 explosion shots (ca. 1.2 kg, 10-15 m depth). A total of 3240 receiver points with 10 Hz geophone groups were all active during the survey. A first processing of the data provided a post-stack time-migrated cube (Stiller et al. 2018;Krawczyk et al. 2019), which was further analysed by the methods described in this paper.

Continuous wavelet transform
Spectral decomposition methods are used to derive spectral amplitudes for given frequency bands or at discrete frequency values of interest. With such techniques, seismic reflection signals can be converted into sets of spectral attributes for the purpose of seismic facies classification. Examples of common spectral decomposition approaches include short-time Fourier transform, wavelet transform, S-transform and combinations of these approaches (e.g. Chakraborty and Okaya 1995;Castagna, Sun and Siegfried 2003;Sinha et al. 2005;Kazemeini et al. 2009;von Hartmann et al. 2012;Wu and Castagna 2017;Huang et al. 2018). In our investigation, we used the continuous wavelet transform (CWT; e.g. Kumar and Foufoula-Georgiou 1997). This method has advantages against the windowed Fourier transform, particularly when spectral information is required with high resolution in both time and frequency as it is the case in the analysis of seismic reflection signals (e.g. Castagna et al. 2003).
A time series signal f (t) is converted into the continuous wavelet transform of the signal W(s, τ ) by use of the integral transform In equation (1), the time series signal f (t) is multiplied by the wavelet function ψ s,τ (t) which, in general form, is defined by Equation (2) implies that a set of wavelet functions ψ s,τ (t) is calculated from a so-called mother wavelet ψ(t) based on variations of the parameters s and τ . Scaling parameter s is used to stretch and compress the wavelet. Parameter τ controls the location of the wavelet along the time axis. Different functions were introduced as mother wavelets, each suited for specific applications. When the sine-like Morlet wavelet (Grossmann and Morlet 1984) is used, the scaling parameter s becomes a period (1/frequency) and equation (1) converts a one-dimensional (1D) time series into a two-dimensional (2D) time-frequency spectrum. A set of complex-valued Morlet wavelets with variable scaling s and shifting τ is determined by Equation (3) represents an oscillating function weighted by a Gaussian envelope. Parameter f 0 defines the basic frequency of the mother wavelet, which relates to the number of oscillations below the Gaussian and, hence, to the length of the wavelet. Prokoph and Barthelmes (1996) introduced the parameter l to describe the length of the mother wavelet: An example for the wavelet transform of seismic reflection data is shown in Fig. 3. The 1D time series (Fig. 3a) is converted into a 2D spectrum as a function of time and frequency (Fig. 3b,c). The 2D sampling in time and frequency is obtained by variable values for parameters s and τ in equations (1) and (4). A shorter and a longer mother wavelet (parameter l in equation (4)) were used in the calculations to generate Fig. 3(b) and (c), respectively. The shorter mother wavelet with l = 5/(2π ) provides higher resolution in time (Fig. 3b), while the longer mother wavelet with l = 10/(2π) provides higher resolution in frequency (Fig. 3c). This effect is related with the uncertainty principle, which said that resolution cannot be arbitrarily small in both time domain and frequency domain (e.g. Prokoph and Barthelmes 1996). In the presented study, we used a relatively short mother wavelet with l = 4/(2π) to emphasize temporal resolution because the analysed reflection signals are relatively short in time.

Self-organizing map
The self-organizing map (SOM; Kohonen 1982Kohonen , 1990Kohonen , 2001 is a neural network type, which is based upon competitive, unsupervised learning. One important field of application is pattern recognition and classification. In the case of seismic facies analysis, seismic signal characteristics are converted into data patterns, which serve as input information for the SOM. After learning, the input data patterns are clustered and the members of each cluster are classified as distinct seismic facies types. Ultimately, the seismic facies types are interpreted as geological facies types. De Matos et al. (2007) suggested a workflow for seismic facies analysis including continuous wavelet transform and subsequent SOM clustering and classification. In this study, we use a modified version of the approach of de Matos et al. (2007) based on the SOM method of Bauer et al. (2008Bauer et al. ( , 2012Bauer et al. ( , 2015. The general workflow includes (a) preparation of data patterns, (b) unsupervised learning, (c) SOM segmentation and (d) classification of all data patterns. Figure 4 illustrates the preparation of data patterns. The seismic reflection horizon of interest is picked, and continuous wavelet transform coefficients are calculated for each trace using a given time window centred at the picked time sample (Fig. 4a). Instead of using all data samples (e.g. de Matos et al. 2007), we apply this data pre-selection step to exclude irrelevant signals and to restrict the information to the specified target horizon.
Data pattern vectors x are generated from the timefrequency domain CWT coefficients w(s, τ ) calculated at single traces within the short-time windows around the picked samples (Fig. 4a,b): The ordering and the meaning of indices for the CWT coefficients are shown in Fig. 4(b). Noteworthy, the choice of the ordering is not critical because of a subsequent normalization procedure: For each component of the data pattern vectors, the mean and standard deviation for all data are determined and used in the normalization of each component (e.g. Bauer et al. 2012). This normalization procedure guarantees that the Type 3 low high CWT absolute value Figure 9 Results of wavelet transformbased signal classification for the EBS horizon. All data patterns which belong to the same seismic facies type are gathered and averaged. inherent information of each data vector component is equally treated within the SOM analysis.
The data pattern vectors x serve as input for the SOM (Fig. 4c). The SOM consists of a two-dimensional arrangement of neurons. Each neuron is associated with a model pattern vector m with the same dimension n as the data pattern vector. Hence, the SOM can be considered as n layers, each representing one component of the model pattern vectors associated with the neurons (Fig. 4c). At the first step of the learning phase, the model vector components of all neurons are initialized by random values (upper row in Fig. 5). During each iteration t of the learning phase, a data pattern vector x j is randomly chosen, and the so-called winning model pattern vector m w with the smallest Euclidean distance (L 2 norm) to the data pattern vector is determined: The model pattern vectors of all neurons i are then updated by with learning rate λ and neighbourhood function η: where r w,i is the distance between winning neuron w and neuron i at the SOM, and parameter σ controls the width of the Gaussian-weighting function. Equation (7) means that all model pattern vectors are modified in a way that they are more similar to the presented data pattern vectors after the operation. This modification is strongest at the winning neuron and its neighbourhood because of the Gaussian-weighting function (equation (8)). With increasing iteration numbers during the learning phase, the random distribution of model pattern vector components is changing into an ordered distribution across the SOM (see three left columns in Fig. 5). This selforganization process is an effect of equations (7) and (8).
Learning continues until changes of the model pattern vectors are getting smaller than a given threshold value.
As a result of the learning phase, similar input data are associated with neighbouring neurons at the trained SOM. We use a total gradient function g to identify regions at the SOM, which have similar model pattern vectors and which are separated from other regions at the SOM by a strong gradient. The total gradient function g i for neuron i is determined by where n is the dimension of the model pattern vector with components m i,k , and x and y represent the axes of the twodimensional SOM, respectively. The development of the gradient function during the learning phase is illustrated in the right column of Fig. 5. The gradient function of the trained SOM (last iteration of the learning phase) is analysed by a segmentation procedure to define clusters of model pattern vectors, which ultimately represent seismic facies types. A watershed segmentation algorithm (Bauer et al. 2012) is used to identify the boundaries between regions with small gradients at the SOM (Fig. 6a,b). For each separated region, a distinct colour code is defined (Fig. 6c). Finally, the gradient function is used to define the colour saturation for each seismic facies type (Fig. 6d). In the application step, each data pattern vector is presented to the trained SOM and the corresponding winning neuron is determined by testing equation (6). The location of the winning neuron at the SOM defines the seismic facies type and related colour code (Fig. 6d).

R E S U L T S A N D D I S C U S S I O N
The above-described workflow was applied for the analysis of the Rotliegend sandstone horizon. The Elbe base sandstone (EBS) is a well-pronounced porous sandstone at the bottom of the Elbe subgroup which is distributed in the marginal areas of the NE German Basin (Gast et al. 1998;McCann 1998). In the E GrSk 3/90 well, the EBS is encountered between 4060 m and 4100 m depth where it shows low, only weakly varying Gamma ray values and decreased P-wave velocity in the sonic log (Fig. 2). The characteristic decrease in P-wave velocity is caused by an increase of porosity (Trautwein and Huenges 2005).
To ensure which part of the waveforms in the processed data represents the EBS layer, we applied finite difference (FD) elastic forward modelling and compared synthetic data with the stacked and time-migrated data at the E GrSk 3/90 well. The modelling was based on a smoothed version of the P-wave velocity from the sonic log (Fig. 7a). Simplified assumptions were made on density (Gardner, Gardner and Gregory 1974) and on Vs by using (Vp/Vs) 2 = 3. The resulting impedance model M1 is shown in Fig. 7(a). We generated another model M2 where the velocity reduction between 4060 m and 4100 m depth associated with the EBS layer was eliminated (Fig. 7a). Synthetic data for the models M1 and M2 were calculated using the FD method of Bohlen (2002). The source wavelet used in these calculations was extracted from the measured data. The centre frequency of the extracted source wavelet was 52 Hz. The differences D of the envelopes between synthetics from M1 and M2 efficiently reveal the effect of the EBS (Fig. 7b) in the time domain data. We used this information to identify the targeted horizon in the stacked and time-migrated data (Fig. 7c).
The horizon was then picked within the entire data cube. We have picked the waveforms at the zero-crossing of the identified reflection signal in correspondence to the modelling results described above. A wavelet transform pattern was determined for each picked signal as illustrated in Fig. 4. One pattern consisted of CWT coefficients for frequencies between 10 and 100 Hz at 2.5 Hz intervals, calculated for each of seven time samples around the travel time pick. The time window was chosen to incorporate variations of the signal spectrum around the picked time. Tests with different time windows provided similar results, indicating a certain robustness regarding this parameter. Altogether, 87,300 wavelet transform patterns were analysed using the self-organizing map (SOM) method of Bauer et al. (2008Bauer et al. ( , 2012Bauer et al. ( , 2015. The CWT patterns were clustered into three different types of reflection signals. Figure 8 shows the distribution of these signal types across the Rotliegend EBS horizon. An important question at this point is, if the results of the signal classification and mapping could be partly influenced by irregularities in the survey geometry, distribution of different source types and parameters used in the processing scheme. Comparison of Figs 1(b) and 8 shows no undesired correlation between the nominal fold and the distribution of signal types along the target horizon. Spectral analysis of signals measured near to the vibroseis and explosion source locations did not show significant and systematic differences in the frequency content. Tests with different processing parameters and different stacking techniques did not show undesired changes in the wavelet transform output and SOM analysis results.
To understand the geological meaning of the different signal types, it is essential to evaluate their defining characteristics. For this purpose, we gathered all individual patterns which belong to the same signal type and determined their representative by arithmetic averaging of each CWT component (Fig. 9). Following this characterization, Type 1 (red colour) shows largest amplitudes around 17-20 Hz, Type 2 (light green colour) around 30-45 Hz and Type 3 (dark green colour) around 35-60 Hz. In general, variations in reflection signal characteristics or seismic facies can be related with a series of lithological properties, including layer thickness, composition of the rock matrix, porosity and pore filling, diagenetic variations, amongst others. In our interpretation, we concentrate on two factors which are known for a strong influence on the frequency content of the reflection signal: thickness variations and variations of seismic attenuation (Figs 10 and 11).
We used FD viscoelastic modelling (Bohlen 2002) based on downhole logging data from the E GrSk 3/90 well to investigate the potential effects of variable thickness and seismic attenuation on the wavelet transform pattern of the Rotliegend EBS reflection signal. For the modelling of thickness variations, we generated a simplified model description with a constant P-wave velocity of 4.8 km/s and a velocity reduction of 0.6 km/s for the EBS. As the measured sonic log shows a more gradual transition instead of a sharp change at the top and bottom of the layer, we used a Gaussian function to describe the low velocity anomaly (Fig. 10a). Because of incomplete density logging data, density was calculated from P-wave velocity (method of Gardner et al. 1974). Seismic attenuation (Q p ) was fixed with an assumed constant value of 80 based on laboratory data for sandstones (Toksöz, Johnston and Timur 1979). Attenuation of S-waves (Q s ) was set equal to Q p in all calculations. Synthetic traces were calculated for models with different thickness (20 m, . . . , 100 m) of the low velocity layer, and CWT patterns were determined from the reflection signal. The results show a strong dependency of the frequency content on thickness (Fig. 10b). Based on this simplified modelling, we can explain the observed CWT patterns (Fig. 9) by predominant thickness of 80-100 m for Type 1 (red colour), around 40 m for Type 2 (light green colour) and around 20 m for Type 3 (dark green colour). In fact, these results are in agreement with a thickness of 39 m at the E GrSk 3/90 well, which is located within facies Type 2, and 22 m at the E GrSk 2/76 well, which is located within facies Type 3 (Fig. 12).
The potential influence of seismic attenuation on the CWT response was based on a model with a 40 m thick zone of reduced P-wave velocity and a variable Qp perturbation around the target layer (Fig. 11a). The modelling results indicate, that a decrease in Qp (increased attenuation) will produce large CWT amplitudes at slightly lower frequencies (Fig. 11b). An increase of Qp (decreased attenuation) leads to a shift of high CWT amplitudes towards slightly higher frequencies. In general, the simplified modelling tests indicate that thickness variations have a much larger influence on the CWT response of the EBS layer compared with reasonable variations of Qp.
Our final interpretation of the seismic facies distribution along the Rotliegend EBS horizon is based on the presented modelling results, borehole information and preexisting knowledge on the geological facies from integration of regional data (e.g. Gast et al. 1998;Gast et al. 2010;Fryberger et al. 2011). Based on this information, we assume a predominantly fluvial facies in the studied area. The redcoloured facies (Type 1) is interpreted as a system of paleochannels (Fig. 12). These parts are related with larger thickness of 80-100 m. If we assume an increase in seismic attenuation (decrease in Q p ), the thickness could be smaller than 80-100 m but still significantly larger than that for the other facies types. The paleo-channels with greater thickness could have been formed within a broader fault zone potentially affected by erosion at the paleo-surface, which would offer a larger space for deposition of sediments. This idea is supported by results of Bauer et al. (2010), where a fault zone was interpreted along those parts of a 2D tomographic model which are crossing the paleo-channels (dashed segment at I-GET profile in Fig. 12). We assume that the paleo-channels developed within a fluvial plain, which shows a predominant thickness of 40 m (Type 2, light green colours in Fig. 12). The fluvial sedimentation was most likely fed from the east, where a topographic swell existed between the NE German Basin and the Polish Trough (e.g. Gast et al. 2010). Thinner sediments (Type 3, dark green colour in Fig. 12) could be partly related with topographic highs and less space for deposition of sediments. Such paleo-topographic highs may have existed as a result of the underlying structure (e.g. Baltrusch and Klarner 1993).

C O N C L U S I O N S
We used the continuous wavelet transform as a full waveform attribute description of reflection signals from a Rotliegend sandstone horizon and evaluated the results in context with geothermal exploration. Self-organizing map (SOM) classification delivered three different seismic facies types which are mapped across the reflector. An important step in our interpretation is the visualization of the underlying average CWT patterns. Characteristic differences were found in the frequency content of each facies type. We used finite difference modelling based on borehole information to investigate the potential influence of variable thickness and seismic attenuation on the CWT response of the analysed horizon. We consider the presented modelling approach as an example for a possible interpretation strategy to improve understanding of seismic facies classification results.
The modelling results show that the observed seismic facies variations can be mainly explained by thickness variations of the sandstone layer. We speculate that thickness variations were controlled by paleo-topographic variations and related variable accumulation space for predominantly fluvial sediments. Thicker paleo-channels were found within a fault zone north of the E GrSk 3/90 and Gt GrSk 4/05 wells. Also, the underlying structure may have influenced paleo-topography and deposition of the sandstones. The results provide constraints for the continuing site development, including the planning of future drilling and geothermal production tests. The lateral distribution of the seismic facies types and their thickness variations can be used to describe the geometry of the sandstone horizon within a new reservoir model. The required assignment of material properties within such a model will be based on our interpretation of paleo-channels. Assumed variations in seismic attenuation together with petrophysical data could further support the definition of material properties. We propose an updated reservoir modelling for Groß Schönebeck based on the new 3D seismic data to investigate the behaviour of the geothermal reservoir during future production activities.
Generally, the presented workflow could be extended if more wells would be available within the seismic study area. In such cases, direct estimation of thickness variations from spectral decomposition and SOM analysis would be possible (Ye et al. 2019).

A C K N O W L E D G E M E N T S
The seismic experiment was financed by the Federal Ministry for Economic Affairs and Energy (BMWi grant 0324065). The authors would like to thank the contractor companies and all people involved in the fieldwork. Support by responsible authorities and local communities is gratefully acknowledged. Ariane Siebert helped to prepare figures. We thank two anonymous reviewers and the Associate Editor for helpful comments and suggestions.