Seismic Lithofacies Distribution Modeling Using the Single Normal Equation Simulation (SNESIM) Algorithm of Multiple-Point Geostatistics in Upper Assam Basin, India

Recently, multiple-point geostatistical simulation gained much attention for its role in spatial reservoir characterization/modeling in geosciences. Accurate lithofacies modeling is a critical step in the characterization of complex geological reservoirs. In this study, multiple-point facies geostatistics based on the SNESIM algorithm integrated with the seismic modeling technique is used as an efficient reservoir modeling approach for lithofacies modeling of the fluvial Tipam formation in the Upper Assam Basin, India. The Tipam formation acts as a potential reservoir rock in the Upper Assam Basin, India. Due to the basin geological complexity and limitation in seismic resolution, many discontinuities in depositional channels in this fluvial depositional environment have been identified using conventional lithofacies mapping. This study combines three techniques to reproduce continuity of the lithofacies for better reservoir modeling. The first is simultaneous prestack inversion for inverting prestack gathers with angle-dependent wavelets into seismic attributes. A cross-plot of P-impedance and VP/VS ratio from well-log data was used to classify the different reservoir lithofacies such as hydrocarbon sand, brine sand, and shale. The second is the Bayesian approach that incorporates probability density functions (PDFs) of non -parametric statistical classification with seismic attributes for converting the seismic attributes into lithofacies volume and the probability volumes of each type lithofacies. The third technique is multiple-point geostatistical simulation (MPS) using the Single Normal equation Simulation (SNESIM) algorithm applied to training images and probability volumes as constraints for a better lithofacies model. These integrated study results proved that MPS could improve reservoir lithofacies characterization.


Introduction
Interpretation of geological shapes and their lithofacies in a reservoir requires a reliable reservoir lithofacies model that is generated with high consistency by integrating the data from various sources (Doyen, 2007). Generally, seismic reflection data and wireline logs are required to identify the reservoir lithofacies (different types of depositional sediments) in many reservoirs.
The lithofacies characterization of a reservoir plays an essential role in identifying prospects essential for development studies.
For many years, lithofacies models have been generated using seismic data and wireline log data by seismic interpretation techniques such as seismic inversion and rock physical studies (Teixeira et al., 2007). The spatial variability of lithofacies in a reservoir can be obtained by integrating wireline logs and seismic data. A widely recognized issue is the ambiguity generally associated with identifying lithology and fluid properties from seismic data (Sengupta and Bachrach, 2007). The main issues associated with seismic interpretation are their limited vertical resolution in geological and petro physical properties (Yu and Ma, 2011). The first source of ambiguity arises from the non-uniqueness of the solutions in seismic inversions in the relationship between the measured seismic amplitudes and elastic properties, as many alternative models can produce the same seismic response. The second source of ambiguity is the link between the elastic properties and petrophysical properties. In rock physics modeling, ambiguity analysis has an essential role in quantitative seismic interpretation (Avseth et al., 2005). However, traditional methods such as seismic inversion and rock physics analysis can understated geological shapes in many petroleum reservoirs, especially in complex reservoirs, leading to ambiguity in lithofacies interpretation.
Statistical modeling methods such as pixel-based two-point simulation and object-based statistical simulation can also create lithofacies models (Deutsch and Wang, 1996;Deutsch and Journel, 1997). However, their reservoir modeling application has been limited due to a lack of conditioning data with wireline logs/seismic data. In two-point geostatistics, lithofacies structures are characterized by using variograms. However, they cannot capture the curvilinear structures and shapes of geological depositions such as meanders, and moreover, the quantity of well-log data is insufficient to obtain a reliable variogram for 3D simulation. On the other hand, objectbased techniques can reproduce lithofacies geometry but limit integrating secondary data and require high CPU usage.
In this study, a reservoir modeling approach is employed that focuses on addressing these challenges using a three-step approach to aid the geoscientist's quest to predict a more accurate lithofacies model by reproducing the curvilinear characteristics (e.g., channels, cross-bedding). It can overcome the limitations of existing geophysical interpretation techniques and statistical methods. This approach involved generating a lithofacies model of the Tipam formation of Upper Assam Basin using wireline logs and 3D prestack seismic data. The methodology of this study integrates the prestack inversion technique and, probabilistic modeling approach with the Single Normal Equation Simulation (SNESIM) algorithm based multiple-point geostatistics method (MPS).
MPS based on the SNESIM algorithm has the benefits of both pixel and object algorithms in lithofacies model generation for complicated geological settings without an explicit non-gaussian model (Hashemi et al., 2014). It is a powerful technique that has proven to be very suitable for simulating the spatial distribution of lithofacies for complex structures (Strebelle, 2002). MPS can be used to model subsurface heterogeneity in complex reservoirs and work as an alternative method to conventional variogram-based approaches that are generally not well suited to simulate complex, curvilinear, continuous, or interconnected structures. The MPS procedure has proved that it is a plausible method for modeling and reproducing hydrogeological underground water channels with the Training image concept (Ti) (Strebelle and Journel, 2002). The concept of multiple-point geostatistics uses the training image (Ti) generated with previous information about facies types and their proportion from the well litho-log and in accordance with the deposition environment. Ti is a strategic concept that is essentially a database of geological patterns/statistics and their associated facies to characterize geological heterogeneity (Fadlelmula et al., 2016). Ti provides multiple geostatistics to the SNESIM algorithm while simulation reproduces the lithofacies channels, which cannot be mapping by conventional modeling techniques. In this technique, no kriging or variogram is involved, but directly receive the probability facies/ patterns from Ti. Another advantage of the SNESIM algorithm is that it reduces the massive demand for memory requirements by avoiding the need to scan the training image for every node simulation, it scans the whole Ti at once as well as and all possible statistics stored as a search tree. This search tree reduces the high CPU usage during the simulation and has no limitations in integrating the secondary data (well-log data and seismic volumes) (Liu et al., 2004).
In this study, we focus on (1) reservoir elastic properties estimating from seismic data, (2) cross plot analysis to identify the lithofacies types, (3) generating probability lithofacies volumes using seismic elastic properties and cross plot analysis through the probabilistic modeling approach, subsequently used as soft data constraints in the SNESIM algorithm, and (4) applying multipoint geostatistics for simulation to reproduce the missing lithofacies and continuity of river deposits that were not mapped in the traditional reservoir modeling of meandering fluvial deposits.

Geological Setting
The Upper Assam Basin is a mature Petrolifeous sedimentary Basin located in the NE Indian subcontinent, covering 57000km 2 . It is part of the Assam-Arakan Basin, and extend up to the Arakan range, Myanmar. The basin limits are defined by the Naga Hills along the south-eastern boundary, the Eastern Himalayas in the north, and the Mikir Hills in the southwest (Sahoo and Gogoi, 2009). The Basin lies between two convergent margins and is categorized by a central NE-SW trending buried "basement ridge," which is a continuation of the exposed Precambrian basement of the Mikir Hills. Paleo-structural analysis indicates this ridge was initiated towards the close of deposition of Barail sediments (Lower Oligocene) and became more pronounced by the end of the deposition of Cenozoic sediments, ranging in age from Paleocene to the Middle-Upper Miocene. Hydrocarbon deposits (oil and gas) are found in almost every stratigraphic horizon. Figure 1 represents the geological location of the study area.
The Upper Assam Basin within the Upper Assam shelf from Miocene to Plio-Pleistocene has many petroleum fields. The important strata with hydrocarbon producing potential are the Tipam sandstone & the Girujan clay formation which belongs to the Miocene age, Barail formation of the Oligocene age, and Sylhet Limestone and Kopili formation from Eocene age. The Tipam and Barail formations are major hydrocarbon producing strata zones in this Basin. The study region is a sandstone deposit within the Upper Assam Basin, India. The Tipam formation is the most crucial lithological unit in this Basin. It is a thick accumulation of fluvial sandstone deposited in a fresh-to brackish water environment in the Miocene age almost throughout the entire Basin and acts as a good reservoir rock for hydrocarbon accumulation in several oil fields Upper Assam. The Tipam formation is subdivided into Upper, Middle, and Lower Tipam. The middle Tipam formation comprises a sand/shale alteration sequence, the Lower Tipam formation is mainly an Arenaceous sequence. The upper Tipam contains an arenaceous sequence (Mallick and Raju, 1995). Figure 2 represents the stratigraphic sequence of the Upper Assam Basin.

Data Set: Wireline Logs and Seismic
Log data from five wells and 3D prestack seismic section were used for the reservoir lithofacies model in the Tipam formation of Upper Assam Basin, India. The field data contain 3D prestack migrated gathers in the offset domain, a total of five wireline logging data, and pre interpreted horizons. The seismic data are zero phase and normal polarity (SEG convention), where an increase in acoustic impedance is displayed as a positive amplitude (black peaks) in the seismic cross-section. The study zone falls within 1600 ms-2300 ms in the time domain (in the seismic section). The seismic offset gathers converted into angle gathers after preconditioning which is required for prestack inversion. Compressional velocity (Vp), Shear velocity (Vs), and Density (Dn), Resistivity logs, Neutron porosity logs are available from five wells. Well#01 and, Well#03 was used in inversion procedure and cross-plot analysis. Well#02, Well#04, and Well#05 were used for quality checks of prestack inversion results and Bayesian classification results. Figure 3 shows raw seismic data in the offset domain (offset gathers), super gathers, and angle domain seismic data (angle gathers).

Methodology
Our methodology was divided into three steps: prestack seismic inversion, probability facies classification, and MPS algorithm. Simultaneous prestack inversion and rock physics analysis were the methods applied to seismic data before applying MPS (Mariethoz and Caers, 2014). Simultaneous prestack inversion was applied to 3D seismic angle gathers to generate seismic attributes . Rock physics analysis via the probabilistic modeling approach is conducted to establish the link between seismic attributes and geological lithofacies identified from the cross-plot analysis (Sen, 2006). Finally, multiple-point geostatistics simulation with the SNESIM algorithm procedure was started with generating the Ti for providing the multiple patterns during the simulation (Mariethoz and Caers, 2014;Mariethoz, 2018). The outcomes of the rock physics analysis (lithofacies probability volumes) used in the MPS approach as soft conditioning and well-logging lithofacies came to the simulation as hard data. This MPS method overcomes the limitation of conventional statistical methods by inferring multivariate data patterns from the training image (Ti). The comprehensive workflow used in the current study is shown in Figure 4.

Estimation of Seismic Elastic Properties from Prestack Seismic Data
As the first stage of the workflow, prestack seismic migrated angle gathers were converted into different geophysical attributes by model-based simultaneous prestack inversions (Russell et al., 2013). Generally, the seismic inversion results are utilized to optimize seismic data and, identify prospects ('sweet spots') in field development studies. Elastic properties are essential for identifying litho-bodies in a reservoir. The seismic data was preconditioned to improve the signalto-noise ratio and converted into an angle domain (angle gathers) for the prestack inversion algorithm (Zhang et al., 2015). Two wavelet extraction methods were used to perform this procedure. First, statistical angle-dependent wavelets were extracted by the autocorrelation method for angle stacks (5°-15°), (15°-25°), and (25°-35°) of seismic prestack angle gathers at two well locations (Gunning and Glinsky, 2006). The well-log reflectivity convolved with these statistical wavelets for generation synthetic seismogram at well locations for correlating with real data. Using the synthetic seismogram, seismic to well-tie operation was conducted at each well location to relate well stratigraphic markers to stratigraphic markers of seismic data and correct the time to depth relation. Later, well-based angle-dependent wavelets were extracted by finding a time-domain operator which shaped the well-log reflectivity that convolved with the real seismic data at every well location (Hampson and Galbraith, 1981;Yi et al., 2013). Initial guess low-frequency models of P-impedance (ZP), S-impedance (ZS), density (Dn), and VP/VS ratio were generated for the seismic data from wireline logs and interpreted horizons. The simultaneous prestack inversion procedure uses these initial models for perturbing iterations with a selected deterministic wavelet to obtain the elastic properties (ZP, ZS, Dn & VP/VS ratio) by achieving a good match between synthetic and real data. The entire seismic inversion approach involves wavelet estimation, well to seismic ties, generating an initial guess model, and seismic inversion based on the reformulated Aki-Richards equation. The prestack inversion was conducted for angle range 5°-35° of seismic data. Figure 5 shows an arbitrary line of ZP & VP/VS ratio from prestack inversion analysis via all five wells. The simultaneous prestack inversion results have been cross verified with well-logging data, including the wells that have

Generating the Probability Lithofacies Volumes Using Logging Data and Seismic Data
The Bayesian approach generates spatial distributions (Probability volumes) by establishing a link between seismic-derived elastic properties and well-derived lithofacies (Avseth, 2000). A cross-plot analysis was conducted using the two wells to define the different lithofacies such as shale (green), water-bearing sand zone (yellow), and prospective-bearing sand represents HC zone (red) with different ranges of P-Impedance (ZP) & VP/VS ratio (Rasaq et al., 2015). The range values of the ZP and VP/VS ratio are shown in Table 1. Lithofacies log for each wellbore location can be estimated from these cross-plot ranges. Figure 6 shows the cross-plot of the three lithofacies and corresponding well litho-logs (Well#01 & Well#03).   After cross-plot analysis, a multivariate Probability Density Function (PDF) was prepared based on non-parametric statistical classification using kernel analysis in the cross-plot space (Loftsgaarden and Quesenberry, 1965). Generally, in rock physics modeling, PDFs are generated using the parametric method. Practical application of the calculation of PDFs using the parametric method is complicated with geophysical data points. It is challenging to assume the parametric family of probability distributions for this kind of data, and precisely estimating the assumed distribution parameters from data is not possible (Ocampo-Duque et al., 2013). A notable disadvantage in the parametric statistical classification is the lack of flexibility. Each assumed distributor's parameters enforce restrictions on the shapes that function may have (Qin et al., 2011). A parametric distribution model may work well for the same kind of geological distribution. A reservoir mostly has ideal characteristics compared with other reservoirs. So, it is challenging to estimate the appropriate parameters in parametric statistical classification. For the kernel analysis, facies proportions for hydrocarbon-bearing sand, water-bearing sand, and shale were estimated from logging data (Węglarczyk, 2018). These proportions of each facies were used in the MPS procedure for generating the Training Image (Ti). Subsequently, these PDFs were used as a priori information in a Bayesian classification to estimate the posterior lithofacies probability from the inverted elastic attributes (Avseth and Mukerji, 2002). Figure 7 shows the PDFs generated from non-parametric kernel analysis. A confusion matrix was used to conduct the QC of different kernels. The confusion matrix measures the mismatch between actual lithofacies (litho-log generated by the cross plot) and estimated Bayesian lithofacies from seismic attributes at well location (Nieto et al., 2013). The diagonal values of the confusion matrix should be large values for the quality of prediction lithofacies and visual inspection on comparison of the actual litho log with the predicted litho log from Bayesian classification. Figure 8 shows the QC measure applied to the Bayesian results using the confusion matrix at well locations. The Bayesian theory is a statistical procedure. Bayesian decision theory provides the primary methodology for classification problem. After proper QC measures of non-parametric PDF at well locations, the Bayesian classification can be applied to the total seismic volume to generate probability volumes for each lithofacies and most lithofacies volume (Loftsgaarden and Quesenberry, 1965;Gramacki, 2018). Using Baye's theorem, the prior knowledge of each lithofacies PDF incorporates seismic attributes to produce the most probability lithofacies volumes (Fukunaga, 1990;Duda et al., 1998). These most probability volumes are used in the MPS procedure as soft data constraints. Figure 9 shows the probability volumes of each lithofacies identified in the cross-plot analysis.

Estimation of Lithofacies Model Using Single Normal Equation Simulation (SNESIM) of Multiple-Point Geostatistics
The MPS procedure involves, two key steps: the first is constructing a searching tree, and the second is the simulation part (Remy et al., 2009). The main steps of the SNESIM algorithm process are as follows (i) Create the Ti with the help of previous geological depositional knowledge with the exact proportion gained from the cross-plot analysis at least 1.5 times the simulation grid. Retrieve the conditional probability distribution from training image P(Z(u)=k|devs(u)).

(viii)
The soft data (probability volumes of each facies) are integrated with conditional probability distribution in simulation.

(ix)
A simulated value is drawn from that conditional probability (Realization).

Generation of Training Image & Construction of a Search Tree
As per the geological knowledge of the study area, Ti is generated with parametric shapes given to each facies using the SGeMS software (Fadlelmula et al., 2016). The geological study suggests that the study area is a part of basin hinge lines with normal faults transacted by reverse faults that preceded the major thrust. The geological interpretation has suggested that the area of study is the vast sedimentation of the Brahmaputra River, having large, braided channels with many meandering channels (Sahoo and Gogoi, 2009). A Ti can be generated with three facies representing shale, water-bearing sand, hydrocarbon-bearing sand based on geological knowledge of the study area, and proportions of lithofacies from the cross-plot analysis. The un-conditional realizations of object-based algorithms were used to generate Ti with sinusoidal channels (sand), elliptical shape geometries (HC), and the remaining area used as background (shale). Both lithofacies can be arranged according to geology depositional direction. Figure 10 shows the generated Ti using parametric shapes.
First, in the scanning of Ti, ti (p) is specified as a value of Ti where p ϵ Gti and Gti is the regular Cartesian grid discretizing the Ti. Generally, ti (p) can be a continuous, categorical, or vectorial variable.
The search template TD consists of a value at central node p and its S nearby nodes Pα (α=1,2,3… S) to capture patterns inside the Ti. Pα is specified as Pα =p+iα (α=1,2, 3…S), where iα is the vector depicting the search template (Tahmasebi and Sahimi, 2016). The search template scans the training image (Ti) to identify the central value patterns (p). the data pattern Pat (p) is a multivariant vector, lithofacies value specified at the central node of the template. Pat (p) = {t(p); t(p + i α ), α = 1,2,3 … S} (1) All these patterns and their possibility categories are stored in the form search tree before simulation. Figure 10. Training image: background represents shale (blue), sinusoidal litho body represents sand (green), and elliptical shape (hydrocarbon).

Generating of Realization (Lithofacies Model) Using the SNESIM Algorithm
In the simulation part, a 3D simulation grid is constructed, and petrophysical data assigned to the closet grid of the wellbore locations. The SNESIM algorithm relies on the idea of simulating each cell value (facies number) or petrophysical property sequentially along a random path, and the simulation of nodes/cells is constrained by cells earlier simulated and hard data (well data) and seismic lithofacies probability volumes (Strebelle, 2002).
The search template TD visits every node of the simulation grid, and the TD is used for retrieving the conditional data event devs (p) at each simulation node u, which is defined as dev s (p) = {z (l) (p + i 1 ) … . . z (l) (p + i s )} Where z (l) (p + is) is an informed nodal value (lithofacies number) in lth SNESIM realization, that value could be either a previously simulated value or an original hard datum value. Any number of uninformed nodal values can be centered at u among the S possible locations of the template TD. Next, the number of replica data patterns are identified from the search tree as data event devs (p). The probability values patterns are retrieved from the search tree, and the value is simulated and assigned a new conditional probability value from it with constraints of soft data (probable volumes) (Hansen et al., 2018). Then the conditional probability function at every node of the simulation grid is calculated for each facies. If the availability of the number of replicates from the search is less than Cmin (the minimum number of replicates), the search is repeated. This procedure works till n ≥ Cmin.
The conditional probability distribution (cpdf) is predicted for a lithofacies value at the template's central node.
F(p: k|dev s (p) = {P(Z(p) = k|dev s (p)) ≈ C k (dev s (p)) C(dev s (p)) where, devs (p)-identified data event with central node value. C(devs(p))-number of statistics similar to devs(p) from the search tree. Ck(devs(p))-number of statistics of center value equal to k.
Integration of seismic probabilities of lithofacies with simulation is explained by (Journel, 2002). Figure 11 shows the 15 th & 22 nd realization of the simulation.

Results and Discussion
In this study, we used a multipoint geostatistical simulation methodology, which allows the incorporation of seismic volumes and wireline data for reproducing the continuity of the lithofacies in the reservoir model of Upper Assam, India. Integrating prior information, such as wireline log data and seismic probability volumes, constrained the simulation, and allowed for a better lithofacies model, and comparison of the results with the Bayesian lithofacies model. The MPS results presented here are an averaged map of the lithofacies produced by averaging the 25 realizations. Figure 12 shows a comparison of a 3D horizontal slice of the lithofacies volume and the MPS lithofacies at 1870-85ms (1815-25m). As per reservoir geology, the Tipam formation is a massive sandstone depositional environment created by rivers processes.  Figure 5 shows that, well#02 and well#03 are 15km apart in this study region, and this clearly indicates the hydrocarbon phase in the Tipam zone at depth range 1815m-1825m (the corresponding slice in Figure 12). The wellbore data (Well-2 & Well-3) clearly shown in this interval represents the hydrocarbon sand phase of the Tipam formation in this reservoir. But the continuity of the zone is missing between these wells. The SNESIM algorithm has reproduced the continuity of lithofacies between the two wells. In the conventional modeling techniques miss the continuity of hydrocarbon sand (left side). The development of hydrocarbon lithofacies by the SNESIM algorithm is shown in Figure 12 (right side). The same improvement in brine-sand can be observed in Figure 12 (red arrows). Figure 13 shows the improvement in lithofacies observed in the lithofacies from Bayesian classification in comparison with MPS simulation at 2050ms-2075ms. This zone is a vast fluvial depositional environment with many meandering channels. As observed in Figure 13, the meandering channel has no proper continuity, but it was developed in MPS model results ( Figure  13 (right)).  The prediction quality of this methodology's confusion matrix is superior that of the confusion matrix of the Bayesian classification. The percentages of the diagonal elements of the confusion matrix have increased due to a reduction in the mismatch between true lithofacies from well log data and lithofacies results from the MPS method. From these confusion matrices, we find that the present methodology warrants reliable lithofacies prediction. The MPS simulation results substantiated that the SNESIM algorithm successfully reproduced the missing lithofacies, which could not be characterized in the conventional Bayesian classification. With real data, the present methodology provides promising alternative to improve lithofacies characterization. Therefore, this study's integrated workflow proved that it could provide a better controlled and reliable lithofacies characterization.

Conclusions
In this work, we integrated seismic modeling techniques with the SNESIM algorithm of MPS. The methodology was successfully applied to prestack seismic data and wireline log data from the Upper Assam. When comparing the MPS results with Bayesian classification results, the SNESIM algorithm results show improvement and continuity in the reservoir's lithofacies. The main stages for developing an updated lithofacies model is involved generating elastic properties from the prestack seismic data using the simultaneous prestack seismic inversion, cross-plot analysis to identify lithofacies types based on P-impedance and VP/VS ratio, conversion of seismic elastic volumes into lithofacies by creating a link between the seismic elastic volumes and lithofacies identification by cross-plotting. The training image (Ti) was generated for the facies identified in cross-plot analysis with geological knowledge of the depositional environment. Finally, the simulation part used the SNESIM algorithm with the assist of multistatistics from Ti and constrain by soft data (lithofacies probabilities) and hard data to reproduce the geological facies. The MPS results were compared with the Bayesian method results. The methodology has been proven an effective technique for improved reservoir characterization. We have combined these approaches to generate a facies model that improves the consistency and continuity of lithofacies in fluvial reservoir scenarios. The quantitative improvement of the MPS results was cross-verified by the confusion matrix.