Application of artificial neural networks for lithofacies determination based on limited well data

Lithofacies de ﬁ nition in the subsurface is an important factor in modeling, regardless of the scale being at reservoir or basin level. In areas with low exploration level, modeling of lithofacies distribution presents a complicated task as very few inputs are available. For this purpose, a case study in the Po ž ega Valley was selected with only one existing well and several seismic sections within an area covering roughly 850 km 2 . For the task of expanding the input data set for lithofacies modeling, neural network analysis was performed that incorporated interpreted lithofacies (sandstone, siltite, marl, and breccia-conglomerate) in a single well and attribute data gathered from a seismic section. Three types of different neural networks were used for the analysis: multilayer perceptron, radial-basis function, and probabilistic neural network. As a result, three lithofacies models were built alongside a seismic section based upon predictions acquired from the neural networks. Three lithofacies were successfully predicted on the section while the breccia-conglomerate was either missing or underpredicted and mostly positioned in a geologically invalid interval. Results obtained by single networks differed from one another, which indicated that a result from a single network should not be treated as representative; thus, the facies distribution for modeling should be acquired from either an ensemble of neural networks or several neural networks. Analysis showed the initial potential of the usability of neural networks and seismic attribute analysis on vintage seismic sections with possible drawbacks of the applications being pointed out.


Introduction
Compared with the other parts of the Pannonian Basin system, the area of the Požega Valley in continental Croatia is a weakly explored region. In this relatively wide structural depression of approximately 850 km 2 , there are only few seismic sections and one exploration well. In order for a complete geologic model of the Neogene-Quaternary infill to be made, lithofacies modeling must be performed, for which a much larger data set is usually needed. For this purpose, artificial neural networks (ANNs) were used, as they had previously proved successful in similar tasks of handling well log and seismic data (Bhatt 2002;Cvetković et al. 2009, Malvić et al. 2010, Bagheripour et al. 2013.

Geologic setting of the exploration area
The Požega Valley is located in the central part of northern Croatia (Fig. 1). It is positioned between highlands that surround it from three sides: Papuk Mt. and Krndija Mt. to the north, Psunj Mt. to the west, as well as Požeška Gora Mt., and Dilj Mt. to the south. Further to the east lies the Đakovo-Vinkovci Plateau, as a slightly uplifted block between the eastern part of the Drava Depression and western part of the Slavonia-Srijem Depression.
Geotectonically, the Požega Valley represents a subdepression of a graben-type subsidence, while a horst-type uplift is associated with surrounding highlands between the Sava and Drava Depressions, according to Najdenovski (1988). Geomorphologically, the subdepression represents a spacious, asymmetrical longitudinal valley, mainly bounded by a 300-m isohypse and covered by Quaternary deposits. The mentioned specific setting allows the formation of an extensive drainage basin, where sediments have been deposited by different mechanisms in alluvial and terrestrial environments, resulting in a lithologically variable sedimentary cover of alluvial sand, gravel, clay, and loess (Ivanišević et al. 2015). The geologic composition of the sediments in this subsided unit is composed of tectonically and stratigraphically different sets of Miocene and Pliocene rocks. There are roughly two rock complexes with different geologic, lithological, petrographic, and genetic characteristics. The first complex contains Paleozoic and Mesozoic rocks that are in the subsurface, usually grouped in one unit, known as pre-Neogene basement (Fig. 2). It includes various metamorphic and magmatic rocks (e.g., gneiss and phyllonite) of Paleozoic age and quartz-chlorite-sericite schists for which Mesozoic age could be presumed (Fig. 2;Pandžić 1979;Pamić et al. 1998). Sporadically, Mesozoic carbonates can be found in the outcrops (Belak et al. 1998), but none were determined in the well Tek-1. The other complex incorporates Neogene and Quaternary rocks, containing sediments and effusives in an age span from Lower Miocene to Holocene. These have traditionally been subdivided into lithostratigraphic units in the rank of formations according to Šimon (1980) and Hernitz (1980). These units are bracketed by E-log marker horizons (Tg, H, G, B, and A; Fig. 2)taken in a broader sense meaning, some of them are actually unconformities, whereas many of them in specific settings can be treated as a chronohorizons (Vrbanac 2002). This conventional correlation is nowadays, obsolete, as there is much to be improved, since very few of these E-log marker horizons follow regional trends (Cvetković 2017). However, it serves well for this (local) example.
Early to Middle Miocene age is represented by rocks of the Vukovar Formation (Fig. 2). This unit, in the surrounding area, is characterized by a diverse lithologyfrom clastic rocks of various grain size, through carbonates, to effusive rocks (Hernitz 1980;Najdenovski 1988;Pavelić 2001;Kovačić et al. 2011;Pavelić et al. 2016). It was deposited during the extensional regime of the forming of the Pannonian Basin with diverse sedimentary environments (Lučić et al. 2001) at the onset of marine transgression during the Middle Badenian (Ćorić et al. 2009). In the well Tek-1, the Vukovar formation is represented by a mixed section of coarse-grained sediments (dominantly breccia) in the basal part and marly limestone in the shallower part of the interval (Fig. 2).
Valpovo Formation is the name given to Lower Pannonian sediments, which were deposited in the calm basin or lacustrine environment (Pavelić 2001). These predominant pelitic sediments are represented by marl and clay-rich limestone. According to Hernitz (1980), at a regional scale, this formation can be missing due to the result of a basin inversion, which occurred during the Sarmatian (Csontos 1995;Horváth 1995;Kováč et al. 1997).
Upper Pannonian and Early Pontian-age sediments correspond to the Vinkovci Formation (Hernitz 1980). They were deposited during the post-rift phase of the evolution of the Pannonian Basin during which the accommodation space was controlled by thermal subsidence (Pavelić 2001;Lučić et al. 2001). The usual thick units consist of numerous layers of marl and mostly fine-grained sandstone deposited in an extensive lacustrine environment (Kovačić et al. 2011).
The Vera Formation relates to approximately Upper Pontian sediments (Hernitz 1980). The lithology is similar to that of the Vinkovci Formation with more coarsegrained (sandstone) intervals, which can be observed in the Tek-1 well (Fig. 2). Geologic column of well Tek-1 modified after Najdenovski (1988) The uppermost part of the Neogene-Quaternary infill belongs to the Vuka Formation of approximately Pliocene, Pleistocene, and Holocene age. Sediments are clay, sand, silt, and gravel with sporadic occurrences of coal. The depositional environment was diverse, ranging from lacustrine (Mandic et al. 2015) to aeolian.
According to Najdenovski (1988), the thickness of the Neogene-Quaternary infill of the Požega subdepression is locally greater than 3,000 m, whereas in the well Tek-1, the total sediment thickness is around 2,000 m. The subdepression represents a closed area of 870 km 2 . Initial exploration did not result in the discovery of any petroleum potential other than small gas occurrences (Najdenovski 1988). Hence, the area remained underexplored with regard to seismic surveys and exploration wellsin total counting, only five seismic sections and one exploration well (Tek-1) in the central part of the subdepression. The interpreted seismic section used for later analysis along with the Tek-1 well is shown in Fig. 3.

Data set analysis and modeling methods
Determination of lithofacies is crucial for the construction of usable subsurface models in any geologic setting. In the case of scarce exploration data, the problem of spatial distribution of lithofacies can possibly be resolved with the help of ANNs. What differentiates ANN from other available methods is the ability to lessen the uncertainty and errors in lithofacies determination through the incorporation of hard data (well logs and seismic attributes) and soft data (predictions) through the training process. The addition of lithology probabilities directly based on well observations will help to reduce uncertainty in reservoir prospecting and qualification (Hami-Eddine et al. 2015) with perspectives to do the same in (semi-) regional studies as suggested in this paper. Input data that can be used for this purpose range from well data and cores to seismic attributes. Seismic attributes are a practical and valuable source of information that is still not always included in regional studies. With only one deep well at our disposal, seismic attributes inevitably gained more significance. They were highly represented in the training process of the ANN, in an attempt to make use of their concealed properties and gain the best possible results.
Today, seismic attributes are a compendium of all the measured, computed, or observed information obtained from seismic data that can greatly contribute to better interpretation of the subsurface. Seismic attributes carry information related to amplitudes, frequencies, and phase. These three are the fundamental attributes of a seismic wave. Multiattribute analysis is already established as a technique to provide researchers with information on seismic anomalies indicating specific geologic structure, different sedimentary environments, fractured zones, and faults. In addition, multiple seismic attributes can be used as direct hydrocarbon indicators (Koson et al. 2014). In this work, numeric values of seismic attributes were extracted and used for further analysis. The advantages and disadvantages of each are described in Table 1. The main constraint of seismic attributes in general is that the ability of good attribute interpretation lies in the quality of the original seismic data. Depending on the direct relation to interpretative significance, seismic attributes are divided into two general classesphysical and geometric attributes. Physical attributes are related to wavepropagation, lithology, and other physical properties. They are furthermore subdivided into pre-and post-stack attributes. Geometric attributes are the way to extract more detailed information regarding the dip, azimuth, and discontinuity (Subrahmanyam and Rao 2008;Pigott et al. 2013). The usage of attributes is mainly focused on 3D seismic cubes; however, they can be computed on 2D sections, as was the case in this study (Fig. 4).

Neural network analysis background
In the past few decades, development in exploration methods and data they provide has led to efforts in trying to find a way to process the large amount of newly available information. A new field of computing and processing has developed for those evergrowing sets of data to be properly analyzed within acceptable time constraints. Thus, the models of artificial intelligence that imitates the principle of the biological neural networks have been proposed as the solution for complex problems. ANNs are a powerful tool that has proved useful in various disciplines, including the geosciences. The basic principle of ANN is so-called template matching, where the Table 1 Representation of theoretical advantages and constraints of 12 seismic attributes used in the paper where bright events signify acoustically strong positive and negative events No.  (2005) (Yegnanarayana 2006). The actual analysis is performed in two steps. First is the iterative process of forming the network, which is called the training process, in which the network is "learning" from the data with the known output. The second phase is the utilization of the network for the purpose of determining the output based on input data or prediction of one or more variables. Every network is represented by its architecture, which differs regarding the number of input and output variables. There is a minimum of three layers in an ANN. The first layer represents the data that is used as input for the analysis. This is called the input layer. The signal travels from input layer to the hidden layer/layers where the processing of the data takes place. Different types of networks are distinguished, based on the connection between the two mentioned layers. The connections between layers can be feedforward and feedback, where every neuron is connected to the previous layer's output, and based on the given result or output, the neuron can be reactivated. The activation is controlled by the activation function. In the feedback model, the weight coefficients are assimilated, thus reducing error, which is why the network is also known as backpropagation network.
In this paper, three different neural network architectures were used to acquire usable lithofacies predictions. Multilayer perceptron (MLP) network is a popular network architecture used in most of the research applications in engineering, mathematical modeling, etc. The MLP network is based on a backpropagation algorithm that has one or more hidden layers and is more successfully applied in classification and prediction problems (Rumelhart et al. 1986).
The radial-basis function (RBF) network is also a commonly used neural network. It is more successfully and frequently applied in solving classification problems than in solving prediction problems (Cvetković et al. 2014). RBF has a hidden layer of radial units (neurons), each modeling a Gaussian response surface. They are also good at modeling non-linear data and can be trained in one stage rather than using an iterative process required by MLP. RBF can learn the given application quickly (Venkatesan and Anitha 2006;Csábrági et al. 2017).
Probabilistic neural networks (PNNs) are useful for automatic pattern recognition, non-linear mapping, and estimation of probability of class membership and likelihood ratios. Regarding data set size, the training speed of PNN can be orders of magnitude faster than the well-known backpropagation algorithm in MLP networks, and still be able to classify unknown patterns with approximately the same success (Specht 1988(Specht , 1990.

Data preparation and ANN training
The first part of the ANN analysis covers the input data preparation process for the training of the network. Lithofacies was interpreted for the entire interval of well Tek-1. The interpretation was based on available well log curves along with cutting and core description from the master log. The analysis itself was carried out in the time domain; thus, the well depths were transformed to two-way time values to be in relation with the seismic data. A simple geologic model was also built in order for the interpreted lithology values to be upscaled. Only four lithofacies categories were used for upscaling as the number of input cases (total of 310) was too small to expect the successful prediction of more lithofacies categories (sandstones, siltite, marl, and breccia-conglomerate). The vertical cell dimension was in correspondence to the sampling of the learning points on the well (e.g., every 4 ms). In total, 12 attributes (Table 1) were extracted from the seismic processing, which represented the input values for the ANN, while the target output was the upscaled lithofacies value.
The networks were trained in StatSoft's analysis software Statistica. Using the "intelligent problem solver" and "custom network designer" options, MLP, RBF, and PNNs were tested alongside ( Table 2). The analysis was made on a thousand networks for each type with a learning rate of 0.01, which resulted in a relatively smooth training graph. Ten of the best networks in the training process were retained alongside one best PNN network. After the review of ANN performance statistics (training and selection error), one MLP network (Table 2) was selected out of the set alongside the PNN network as standalone networks for predictions. An ensemble (ENS) of 10 networks (excluding the PNN) was also built so the networks could compete among themselves for determining the output. Upscaled lithology values predicted from MLP, PNN, and network ENS can be observed in Fig. 5.

Results and discussion
Trained ANNs were used to predict the lithofacies values across one of the seismic sections within the Požega subdepression. Prediction was performed on every 10th seismic trace in the section with a depth resolution of 4 ms. The ANN methods gave different results. The first one was derived from the prediction of the standalone MLP ANN (Fig. 6a), the second from PNN ANN (Fig. 6b), and the third from an ANN ENS of 10 different networks (Fig. 6c). The MLP network, although it showed the best results in the training process, predicted a large spread of sandstone facies, which is not in correspondence with the general geologic setting that can be observed in the outcrops. Even though the fitting of the interpreted upscaled versus predicted facies is much more accurate for sandstone than in other ANN analyses (Fig. 5), such a high sandstone content predicted by MLP network should not directly be treated as an erroneous result but rather as one of the probable results. The predominant facies that was predicted with the PNN was marl, which could be expected regarding the training result observed in Fig. 5, with underestimated sandstone lithofacies. Coarser-grained facies were determined in the Vukovar and in the Vuka Formations, which can be related to the regional geologic settings described in Pavelić (2001) and Kovačić et al. (2011). The ENS of ANNs presented similar results as the PNN, with slightly less Central European Geology 60, 2017 sandstone predicted. What was evident in all three cases is the inability of the ANNs to predict the breccia-conglomerate lithofacies within the stratigraphic interval in which they are expected (e.g., Vukovar Formation). Predictions of this lithofacies are either sporadic and not geologically valid in terms of regional settings (Fig. 6a, cbrecciaconglomerates determined in formations of approximate Late Miocene age) or almost absent (Fig. 6b). From the mathematical point of view, the prediction of the brecciaconglomerate lithofacies in stratigraphically erroneous positions is valid, as a combination of seismic attributes could occur in such intervals. Regarding this, the learning process should possibly include more variables, not only seismic, which could be spatially defined across the seismic section.
To improve the obtained results, a few steps could possibly be taken. Results such as those regarding the breccia-conglomerate lithofacies are a result of a small number of input cases through which the network was learning the properties of the lithofacies. Increasing the number of input cases from wells outside the exploration area, but from similar geologic environments (e.g., wells from the Sava, Eastern Drava, and Slavonia-Srijem Depressions), would be valuable. Prediction could be further enhanced by training the neural networks for each stratigraphic interval separately, which would exclude the possibility of predicting breccia-conglomerates in the stratigraphically inappropriate interval. This would only be possible if the learning data set were to be increased, by the number of cases per stratigraphic interval in the current data set (e.g., only 14 cases within the Valpovo Formation, which is visible on Fig. 5). The breccia-conglomerate problem could also be resolved by adding an additional categorical variable, which would determine the stratigraphic properties of a sample point, i.e., defining the facies constricted to the zone of Middle Badenian based on well picks and subsequently by horizon interpretation. In this way, a neural network could possibly recognize that the breccia-conglomerates occur only in the Vukovar Formation and not predict them in other formations, regardless of attribute values.
The large difference of predicted lithofacies values across the section in standalone ANNs suggests that one single solution should not be regarded as a representative one.

Conclusions
The usage of ANNs has shown that it is possible to successfully train them and acquire predictions based on attribute values across the seismic section. Although the presented results differ from each other, they show a certain relation to the regional geologic settings. A result from a single neural network (standalone) should not be treated as representative, as the analysis showed that predictions from different networks vary. The breccia-conglomerate facies proved to be the most underpredicted, as it is either absent or predicted in a small number of cases that are not at all geologically valid.
The number of lithofacies that can be predicted in this way varies in terms of the number of available cases; the presumption is that more cases could certainly yield a more successful prediction of lithofacies. Further testing of such an application should be done in several directions. First is the expansion of the input data set with well and seismic data from geologically similar environments, e.g., the Sava, Eastern Drava, and Slavonia-Srijem Depressions. Second, additional variables should be added to restrict the prediction of certain facies within a stratigraphic interval, so that the resulted prediction is more geologically valid. As the prediction from one network could not possibly be representative, results from several networks should be obtained for the building of a geologic model.