Integration of Adaptive Neuro-Fuzzy Inference System, Neural Networks and Geostatistical Methods for Fracture Density Modeling

Image logs provide useful information for fracture study in naturally fractured reservoir. Fracture dip, azimuth, aperture and fracture density can be obtained from image logs and have great importance in naturally fractured reservoir characterization. Imaging all fractured parts of hydrocarbon reservoirs and interpreting the results is expensive and time consuming. In this study, an improved method to make a quantitative correlation between fracture densities obtained from image logs and conventional well log data by integration of different artificial intelligence systems was proposed. The proposed method combines the results of Adaptive Neuro-Fuzzy Inference System (ANFIS) and Neural Networks (NN) algorithms for overall estimation of fracture density from conventional well log data. A simple averaging method was used to obtain a better result by combining results of ANFIS and NN. The algorithm applied on other wells of the field to obtain fracture density. In order to model the fracture density in the reservoir, we used variography and sequential simulation algorithms like Sequential Indicator Simulation (SIS) and Truncated Gaussian Simulation (TGS). The overall algorithm applied to Asmari reservoir one of the SW Iranian oil fields. Histogram analysis applied to control the quality of the obtained models. Results of this study show that for higher number of fracture facies the TGS algorithm works better than SIS but in small number of fracture facies both algorithms provide approximately same results.


INTRODUCTION
The word "fracture" is used as a collective term representing any of a series of discontinuous fractures in rocks such as joints, faults, fractures and/or bedding planes. In naturally fractured reservoir fractures have a significant effect on hydraulic properties of reservoir. When fractures are open, they act as pathways for hydrocarbon production and may even transform a very low permeability reservoir into highly productive zones. When cemented they act as barriers to hydrocarbon flow hindering the motion of hydrocarbon toward the wells (Haller and Porturas, 1998;Khoshbakht et al., 2009). Therefore in modeling fractured reservoirs, understanding fracture properties is very important. Natural fractures can be identified and evaluated by several techniques, with the most common being core analysis, well log analysis and pressure transient testing. Since mid-1980s and introduction of image log tools, the process of fracture detection and characterization of fracture properties; such as dip, dip direction and fracture density; has become less problematic (Serra, 1989;Tokhmchi et al., 2010). As fracture modeling with an inadequate volume of data can lead to misleading interpretations any direct or indirect techniques which increase the knowledge of fracture properties is highly valuable (Tokhmchi et al., 2010).
In this work, we tried to find the best method to predict fracture density from petrophysical log data. NN and ANFIS methods were used for this purpose. ANFIS model was previously done by Ja'fari et al. (2012) and the results are available. A same dataset was used for NN to predict the fracture density. A simple averaging method was also used to average both model outputs. The results show that output of this average method has the best correlation with real fracture densities. So this method was used to predict the fracture density in all other wells of the field. Finally the fracture density was modeled in the inter-well region using Sequential Indicator Simulation (SIS) and Truncated Gaussian Simulation (TGS) and histograms were used to quality control of models.

METHODOLOGY
In this work, geostatistical and artificial intelligence methods were integrated to create a fracture density model in entire of the field. A brief description and background of each method was described here.

Adaptive Neuro-Fuzzy Inference System
Neuro-fuzzy modeling is a technique for describing the behavior of a system using fuzzy inference rules within a NN structure. Using a given input/output data set, adaptive Neuro-Fuzzy Inference System (ANFIS) constructs a FIS whose MF parameters are tuned using a back propagation algorithm (Matlab User's Guide, 2007;Labani et al., 2010). So, the FIS could learn from the training data. FL and ANFIS were used by different authors to predict target parameters from a set of input data. Gokceoglu et al. (2004) used neurofuzzy model for modulus of deformation of jointed rock masses. Kadkhodaie-Ilkhchi et al. (2009) used a committee fuzzy inference system to predict petrophysical data from seismic attributes. Labani et al. (2010) used a committee machine with intelligent system to predict NMR log parameter from petrophysical log data.

Neural Network
NN is an intelligent tool for solving non-linear problems. Back propagation, or propagation of error, is a common method of training Artificial Neural Networks (ANN) to learn how to perform a given task. It's a supervised learning method; it means that it requires a set of training data that has the desired output for any given input. The network computes the difference between the calculated output and corresponding desired output from the training data set. The error is then propagated backward through the net and the weights are adjusted during a number of iterations, named epochs. The training stops when the calculated output values best approximate the desired values (Bhatt and Helle, 2002;Labani et al., 2010). NN are used widely in petroleum geoscience to predict different target parameters and also to model a parameter all over a field. FitzGerald (1999) used ANN to predict fracture frequency from petrophysical logs. Ouenes (1999) introduced a new approach in fractured reservoir characterization which uses fuzzy logic and NN. El Ouahed et al. (2005) developed a 2D fracture intensity map obtain a better result by combining results of ANFIS and NN. The algorithm applied on other wells of the field to obtain fracture density. In order to model the fracture density in the reservoir, we used variography and sequential simulation algorithms like Sequential Indicator Simulation (SIS) and Truncated Gaussian Simulation (TGS). The overall algorithm applied to Asmari reservoir one of the SW Iranian oil fields. Histogram analysis applied to control the quality of the obtained models. Results of this study show that for higher number of fracture facies the TGS algorithm works better than SIS but in small number of fracture facies both algorithms provide approximately same results. and fracture network map in a large block of Hassi Messaoud field, using artificial NN and fuzzy logic. Darabi et al. (2010) showed the applicability of ANN and fuzzy logic in characterizing Parsi naturally fractured reservoir.

Variograms
When modeling a reservoir with a pixel-based technique, one has to resort to semivariograms to model the sizes and spatial distributions of the parameter. The user usually is faced to a choice of standard (e.g., spherical, exponential, Gaussian) and nonstandard (e.g., fractal, linear, sinusoidal) variograms. The idea is to model existing data (e.g., core, log) with such variograms until the ''best model-fit'' provides the variogram structure (often nested of several components) that is used for reservoir modeling. Because data are often complex and the sample semivariogram can rarely be matched accurately, the problem arises on which semivariogram model should be used (Seifert and Jensen, 1999). In this study, spherical, semivariograms were used for our works.
There are many hydrological and geological processes where heuristic considerations suggest that different values of the variable (low and high) possess different variograms. One way of capturing different variograms for high and low values is to use indicator techniques (Journel, 1983(Journel, , 1993. Indicator semivariograms are similar to traditional semivariograms, except that they are calculated on indicator values rather than the actual value of the variable of interest. Indicator values indicate whether the value of the variable is above (indicator value = 1) or below (indicator value = 0) a particular threshold. The threshold is usually given as the percentile of the univariate distribution of the variable. Since indicator semivariograms can be calculated for a number of different thresholds, they allow for a different spatial structure (i.e. semivariogram) at each threshold (Western et al., 1998).
Different authors provided a comprehensive literature on geostatistical methods and their application in petroleum engineering, like Deutsch (2002). Kelkar (2000). Yarus and Chambers (2006) provided an overview for application of geostatistics for reservoir characterization. Gringarten and Deutsch (1999) introduced a methodology for variogram interpretation and modeling for improved reservoir characterization. Liu and Journel (2007) presented a package for geostatistical integration of coarse and fine scale data. Deutsch (2006) illustrated a SIS program for categorical variables with point and block data. Gringarten (1998) offered a computer code for stochastic simulation of fractures in layered systems.

Sequential Indicator Simulation (SIS)
The SIS method can be used for the stochastic modeling of discrete (e.g., rock types) and continuous attributes (e.g., porosity and permeability). In this method, each attribute to be modeled is described through a binary indicator variable that takes the value ''1'' if that attribute is encountered at a given location, ''0'' if not. Each indicator variable is in turn defined by its average frequency and a semivariogram that characterizes the spatial continuity. Following a random path through the three-dimensional grid, individual grid-nodes are simulated, one after another, using constantly updated, thus increasing size and conditioning datasets. The conditioning includes the original data (e.g., well data) and all previously simulated values within a specified neighborhood. This ensures that closely spaced values have the correct short scale correlation. In other words, in this simulation approach, a grid-node is selected randomly and simulated with reference to the original conditioning dataset (e.g., well data). In the next step, another grid-node is selected randomly, and the variable is simulated using the newly generated CCDF (Conditional Cumulative Distribution Function), which is now increased in size by one value. In this way, each node is simulated (Deutsch and Journel, 1992;Seifert and Jensen, 1999).

Truncated Gaussian Simulation (TGS)
Truncated Gaussian Simulation method has been first designed to provide stochastic images of sedimentary geology, mostly in fluvio-deltaic environments. The basic principle of the method is to replace the handling of the geological description, poorly designed for calculations, by the handling of a Random Function with a multigaussian distribution, for which geostatistical simulations are used routinely. The rock type provided by the simulation depends on the value provided for the Gaussian random function. After definition of thresholds, the rock type is chosen depending between which threshold is the value of the Gaussian random function. The thresholds are given by the proportions of each rock type. These proportions will vary in space, while the properties of the Gaussian random function stay the same (stationary). This method (Matheron et al., 1987) relies on the truncation of a single Gaussian Random Field (GRF) in order to generate realizations of lithofacies. The main feature is the reproduction of the indicator variograms associated with the lithofacies and the hierarchical contact relationship among them. This method is adequate for deposits where the lithofacies exhibit a hierarchical spatial distribution, such as depositional environments or sedimentary formations. The procedure to obtain lithofacies realizations using TGS is described as follows: -establish the lithofacies proportions and their contact relationships; -using the truncation rule, perform variography of the lithofacies indicators through the determination of the covariance function of the underlying GRF; -simulate the GRF at the data locations conditionally to the lithofacies coding. This step is performed using the Gibbs sampler algorithm (Geman and Geman, 1984); -simulate the GRF at the target locations using the values generated at the previous step as conditioning data; -truncate the realizations according to the truncation rule.
The TGS here used as an alternative method for fracture facies modeling and the categorical variograms of SIS are used for this method.

FIELD DESCRIPTION
This field located in SW of Iran has three major faults. 15 wells were drilled in this field until now. For 12 wells the petrophysical log data are available and 2 wells have image logs. 4 productive zones in fractured Asmari formation exist in this field. Variography and data analysis were done separately for each fracture facies. First, the structural model was built and then fracture modeling was done in the framework of this structural model. Alavi (2004) collects all lithostratigraphic units of the Zagros fold-thrust belt data and published the perfect description of formations. Asmari (Oligocene to lowermost Miocene), medium-bedded to thick-bedded, locally shelly or oolitic, nummulites-bearing limestones (grainstone, packstone, wackestone) shoaling upward above a thin basal conglomerate from fine-grained (low-energy) deep-marine marly limestone to high-energy shallow-marine skeletal grainstone; composed of a number of sequences; a unconformity-bounded, highly prolific reservoir; interpreted as transgressive-regressive foredeep facies of the proforeland basin (Khoshbakht et al., 2009).

ESTIMATING FRACTURE DENSITY
The relationship between fracture densities and well log data including sonic log (DT), deep resistivity log (Rd), neutron porosity log (NPHI) and bulk density log (RHOB) are shown in the crossplots of Figure 1. The data used in Figure 1 are normalized between 0 and 1. As shown, a strong and direct relationship between fracture density and bulk density log is seen (CC = 0.5855). Sonic log shows a strong and inverse relationship with fracture density. The dense rocks have a high potential for fracturing. The increase in density of a rock causes the sonic log to decrease. The neutron and resistivity logs show a weaker relationship with fracture density.
Data scaling is necessary here for two reasons. First, it is desired to account for essential variability in the filtered log data and without some type of scaling process, those logs with the largest original variance would dominate the subsequent analysis. Second, it is desired to have all logs measured in similar units because it will be easier to compare them in the next models and the analysis will not be biased towards those with higher absolute values. In this study, a linear scaling method that maps the maximum log value to 1 and the minimum log value to 0 was used. The linear scaling has the following form: (1) where z i is the scaled value, x i is the original value, x min is the minimum log value and x max is the maximum log value. Crossplots showing relationship between fracture density and a) sonic, b) bulk density, c) neutron porosity, d) logarithm of deep resistivity.

NF Model
The formulation between input and output data is performed through a set of fuzzy if-then rules. Normally, fuzzy rules are extracted through a fuzzy clustering process. Subtractive clustering (Chiu, 1994) is one of the effective methods for constructing a fuzzy model. The effectiveness of a fuzzy model is related to the search for the optimal clustering radius, which is a controlling parameter for determining the number of fuzzy if-then rules. Fewer clusters might not cover the entire domain, and more clusters (result in more rules) can complicate the system behavior and may lead to lower performance . Depending on the case study, it is necessary to optimize this parameter for construction of fuzzy model. The clustering radius of the system is changed to find the best clustering radius for the system. After selecting the best clustering radius the data are predicted using this radius. Clustering radius of 0.3 was selected as the best radius for our data prediction with low error and high correlation and also very rapid regard to clustering radius of 0.2 and 0.1. The mean square error of the fuzzy model was 3.7 x 10 -4 and Figure 2 shows a cross plot of NF predicted and real fracture densities.

NN Model
A simple feed forward, back propagation error algorithm was used to construct the NN model. A same dataset as ANFIS model was used for the NN model. The learning function of the network was Levenberg-Marquardt algorithm (TrainLM). A same dataset as NF model was used for NN model. Constructed network had a hidden layer with 30 neurons which used Tansig as transfer function and the output layer with one neuron which used Purelin as transfer function. The performance mean square error of the model was 7.82 x 10 -4 for fracture density. Figure 3 shows the correlation between predicted and real data in the training dataset. This high correlation shows the ability of the network to find the relationship between input and output data. We tried to find the best method for fracture density estimation, so we made a simple average of both models outputs and plot the average results versus real data, too. Figure 4 shows the result of this cross plot. R 2 of 0.9879 shows this is the best fit model among the above models. Crossplot of real and estimated average of fracture density.

Figure 5
Real and predicted fracture density versus depth. Average fracture density map for reservoir interval.
So we used this method to estimate fracture density in other wells of the field. In Figure 5, the value of real and predicted fracture densities are plotted versus depth. This plot shows the high similarity of predicted and real fracture density versus depth and also the great ability of the proposed method to predict fracture density. The mean square error for this method is 2.8 x 10 -4 .
The result of image logs shows that the number of fractures per meter of length of boreholes changes from zero to eight fractures per meter. This can be seen in the predicted fracture density data, too. Therefore 9 facies for fractures were made and each facies indicates a fracture density. For example, facies zero consists of zero fracture per meter, facies one consists of one fracture per meter and so on. These facies are called fracture facies. Geostatistical simulation variograms were calculated and modeled for each fracture facies, separately. The variograms were calculated in a major direction (where the sample points have the strongest correlation), a minor direction (perpendicular to the major direction) and vertical directions using Petrel software. Variography is composed of three steps: determination of experimental variogram, fitting of a model to this variogram and determination of variogram parameters. At first deferent point pairs with an identical separation distance (lag distance) in an identical azimuth were determined. Then, the mean square errors at this points (or variance according to each lag distance) were calculated. This method is repeated for each lag distance. Finally to obtain an experimental variogram, these variances were plotted versus lags. In order to extract the variogram parameters, a model should be fitted to this experimental variogram. Different existing models were tested and finally the model that best fitted the experimental variogram is selected. Before calculating the variograms, the fracture facies proportions in each subzone is determined. Figure 7 demonstrates the proportion of fracture Fracture facies proportion in subzones of 1-4, 2-5, 5-4 and 6-5 are shown in a-d respectively.

Figure 8
Variogram of two fracture facies in major direction. facies in the original data of 4 subzones. Figures 8 and 9 show 2 calculated variograms by using spherical modeling. Gray points on the variogram show the semi-variance calculated for each lag distance. There is a histogram in back ground of each variogram. The histograms indicate the number of "point-pairs" used in the variogram calculation. It shows that each point plotted in the variogram is represented by x number of point pairs. The small histograms bars towards a lag distance illustrate that the variance calculated for that lag distance is not represented by enough point pairs. The variograms have two curves, the gray one which is the best fit for the data and the blue one used for modeling is standardized to sill of one. The parameters of variograms are written below them.

Figure 10
Fracture density model using SIS algorithm. a) Cross section, b) longitudinal section of the model. of stochastic simulation that creates multiple equiprobable models. However, the probability of each model is equal to the other one but the characteristics of them are different. Histograms can be used to control the quality of resulted models (Petrel User's Guide, 2009). Figures 10 and 11 show the results of fracture density modeling, using SIS and TGS algorithms. Figures 12 and 13 show the histograms of upscaled well logs and modeled properties for fracture density. Similar distribution of histograms for each model would verify the accuracy of obtained models.   Histogram of upscaled and modeled fracture density using SIS algorithm.  Histogram of upscaled and modeled fracture density using TGS algorithm.   The results show that histograms of upscaled and modeled fracture density are more similar in TGS algorithm than in SIS algorithm. So in these proposed histograms obtained with TGS work better than the one obtained by SIS.
In another attempt to model fracture densities, the number of fracture facies is reduced by classifying the fracture densities into three classes as class one, two and three containing 0-2 fractures per meter, 3-5 fractures per meter and more than 6 fractures per meter, respectively. Variography analysis and geostatistical modeling of fracture densities were done, using SIS and TGS algorithms. The results of modeling and the histograms of models are shown in Figures 14 to 17. The histograms of modeled and upscaled fracture densities are more similar than the previous modeling results with nine fracture facies. This point verifies that decreasing the number of fracture facies may lead to increase the accuracy of modeling results, when using SIS and TGS algorithms.

CONCLUSIONS
In this study, artificial intelligence tools and geostatistical methods were used to model fracture density in subjected reservoir. The results show that integration of Neural Networks, ANFIS and geostatistical methods could provide very useful data for reservoir characterization and development. The main conclusions of this paper are: -ANFIS and Neural Networks have a great ability in determining relationships between a series of petrophysical log data and fracture density as target; -similar results of ANFIS and Neural Networks verify the trueness of the fracture density prediction; -simple averages of both ANFIS and Neural Networks lead to a better correlation between real and predicted data, so this method can be used for predicting fracture density in this reservoir; -integration of ANFIS and Neural Networks gives improved results for this dataset and provide the required data for the next step of fracture density modeling; -the numbers of point-pairs included in variograms are sufficient to obtain reliable results. The nugget effects in all variograms are reasonable. So, we conclude that variograms can determine the spacial variation of fracture data; -histograms of models show the preference of TGS algorithm to SIS algorithm for this dataset. Even if both models are different in details, histograms show a similar distribution for upscaled and modeled property for both algorithms. This similarity in results verifies the trueness of modeling; -histograms show that a decrease in the number of fracture facies lead to an increase in the accuracy of obtained models.  Fracture facies Figure 16 Histogram of upscaled and modeled fracture density using SIS algorithm.

Figure 17
Histogram of upscaled and modeled fracture density using TGS algorithm.