Pedram Masoudi
email: masoudip@alumni.ut.ac.ir

Behzad Tokhmechi

Alireza Bashari

Majid Ansari Jafari

Bandar Abbas

Embayment Hinterland Dezful

Zagros F High

Raz Ak

Identifying productive zones of Sarvak formation by integrating outputs of different classification methods

Keywords: Sarvak formation, carbonate reservoir, productivity, pay zone, zoning, data fusion Q1

Sarvak formation is the second major carbonate reservoir in Iran. There are several geological, petrophysical and geophysical investigations which have been carried out on this important reservoir. In this work, Sarvak is studied to find productive zones. At first, four different methods were used to identify producing intervals from well log data and well test results. Then, final zoning is generated by integrating outputs of these four methods. One of them is the conventional cut-off-based method; the other three methods are based on flow equation, Bayesian and fuzzy theories. Thereafter, by considering the classification correctness rate of each classifier in each well and technique of majority voting, a unique zoning for Sarvak formation is presented. Based on the final zoning, the whole Sarvak interval is divided into seven zones. Three of them are classified as oil producing zones, two of them cannot be classified as conventionally producing zones, and the remaining two are water producing. Zone number 2 not only has the highest production rate, but also is the most homogeneous zone among the productive zones. The novelty of this research is using well test results in defining productive classes, which improves the certainty of classification in comparison with previous works that were based on core analysis and log data.

Introduction

Sarvak is a carbonate reservoir rock, deposited during Middle Cretaceous. Probably the first geologists who published a paper about Sarvak formation were James and Wynd. In this paper, two members of Sarvak formation were identified [START_REF] James | Stratigraphic nomenclature of Iranian Oil Consortium Agreement Area[END_REF]. There are lots of investigations on Sarvak formation (and its equivalents, especially Mishrif) from different points of view. Here, some of the recent research works, most relevant to this work, are presented.

In 2007, by integrating seismic and well log data, a 3D static model of the Mishrif reservoir (upper Sarvak) was created in Sirri district, the Persian Gulf. According to the static model, it was concluded that a gradual transition exists between Khatiyeh (a basinal deposition facies, equivalent to Q3 middle Sarvak) and Mishrif facies but the transition zone varies in thickness (Bashari 2007b). Bashari has completed another Q4 investigation about the Mishrif reservoir in the Reshadat oil field, the Persian Gulf. In this study, seismic, petrophysical and petrographical data were combined to create a static model of Mishrif, and then, reservoir quality, thickness variation and lateral facies changes of Mishrif were discussed (Bashari 2007c).

Q5

Another valuable work was carried out in the same year by Taghavi et al. As an output of this paper, upper Sarvak is classified into eight different flow units in the Dehluran field, SW Iran. This classification is obtained from evaluating the heterogeneity of petrophysical parameters, after assessing the heterogeneity and distribution of geological variables. Out of these eight distinguished flow units, two (FU3 and FU6) belong to grain-supported shoal/reef deposits with good reservoir qualities. FU1 and FU4 are related to leached lagoonal deposits, the main pore types of which are vuggy or micropores, and represent medium reservoir quality, and the remaining four are non-producing flow units (FU2, FU5, FU7 and FU8) (Taghavi et al 2007).

In 2009, a comprehensive study was conducted by Shahvar et al about predicting permeability in Sarvak formation from different aspects. They predicted a flow zone index (FZI) using fuzzy logic and well log data. Then, the FZI approach was applied for reservoir rock typing and defining hydraulic units. Thereafter, an artificial neural network (ANN) was utilized to estimate the permeability of Sarvak formation in each of its hydraulic unit separately. Permeability was also identified by a conventional regression-based method. Finally, they concluded that the conventional method is not as robust as the intelligent one in predicting the behaviour of heterogeneous carbonate reservoirs (Shahvar et al 2009).

There are some other research works on this reservoir rock which have studied several properties of Sarvak regarding several viewpoints. In 2009, Al-Ameri et al analysed the biomarkers, palynofacies and pyrolysis results of the Mishrif reservoir in Zubair and both Rumaila North and South fields, Kuwait area, and found the source rock of this reservoir in Upper Jurassic Formations [START_REF] Khafaji | Petroleum system analysis of the Mishrif reservoir in the Ratawi, Zubair, North and South Rumaila oil fields[END_REF] Here, a zoning on the whole Sarvak interval is presented regarding the production rate. The article, published by Taghavi et al (2007), is closely related to the current work. They have classified only upper Sarvak into different flow units by analysing geological and petrophysical raw data, whereas here, the whole Sarvak is classified into different producing zones due to production rates, derived from well tests.

Geological setting

This investigation is applied on the Sarvak interval of six wells, drilled on an anonymous field in the Abadan Plain, SW Iran (figure 1). The gross interval of Sarvak in these six wells varies approximately from 590 to 670 m, and is restricted to the Khazhdumi and Laffarn Formations from bottom and top, respectively (figure 2).

Sarvak formation (Upper Albian to Upper Turonian) is an important carbonate reservoir of Bangestan Group, and stratigraphically is equivalent to Mauddud, Khatiyah and Mishrif sequences in the Persian Gulf (Bashari 2007a). Simplified stratigraphy of Cretaceous is presented in figure 2. A detailed stratigraphic study has revealed that Sarvak is deposited on a shallow carbonate shelf setting. Moreover, generally, younger parts belong to shallower environments Bashari (2007a).

in comparison to older parts. While Sarvak is overlaid on the Kazhdumi Formation with a transitional contact, its upper boundary is an unconformity, interpreted as uplift [START_REF] Ghabeishavi | Microfacies and depositional environment of the Cenomanian of the Bangestan anticline[END_REF].

Methods and dataset

The dataset of each well includes well log data, core tests and production test results. Due to the inconsistency of well logs, those logs which are available in all the wells are selected and the others are ignored. The whole dataset and utilized well logs are introduced in table 1.

In this research, productive zones of the Sarvak interval are determined five times by four different methodologies. The first one is a conventional net pay determination procedure, which utilizes cut-off values on petrophysical data. The second approach is based on the flow equation. The third method is a newly developed procedure, based on the concept of conditional probability (Bayes theory). And the last methodology is based on fuzzy theory.

Conventional net pay determination

Conventionally, net pays are determined by applying cutoff values on porosity, shale percentage and water saturation values. There are a few procedures to find the cut-off values of petrophysical data. For a comprehensive study about the conventional methodology, respected readers are referred to Worthington and Cosentino (2005), [START_REF] Jensen | Some statistical issues in selecting porosity cutoffs for estimating net pay[END_REF] and Worthington (2010). Here, cut-off values provided by National Iranian Oil Company (NIOC) reports are used to Q8 determine net pay zones through the Sarvak interval. The output of this methodology is crisp, i.e. productive or nonproductive.

Diffusivity equation

A newly published procedure to determine productive zones is based on the flow equation. In this method, the diffusivity equation is solved by the Ei-function solution, and then rearranged into the form of relation (1). The index of productivity in this relation is the ratio of the flow rate to the pressure differential, which can be easily achieved by reversing the left-hand side of relation (1). For a detailed study about this methodology, the respected readers are referred to Masoudi et al (2011).

p(r, t ) Q 0 = 70.6 μ 0 kh × in 948 ϕμ 0 c t r 2 kt , ( 1
)
where p(r, t) is the pressure at radius r from the well after t seconds in psi, t is the time (h), k is the permeability (mD), Q 0 is the flow rate (rb d -1 ) and h is the pay zone thickness (ft). Using the Ei-function solution in solving the diffusivity equation incorporates approximation error. To keep this error less than 0.25% input of logarithm should not exceed 0.01 [START_REF] Ahmed | Reservoir Engineering Handbook[END_REF].

Bayesian classifier

This methodology is based on Bayes theory. The algorithm Q9 used here is a parametric Bayesian classifier, the discrimination function of which is introduced by [START_REF] Duda | Pattern Classification[END_REF].

• One of the wells which contains all the probable conditions is utilized as a training well to train the Bayesian classifier through it. The trained Bayesian classifier will be applied on all other wells too.

P(npv) = n i=1 P(npv|d i ) × P(d i ), ( 2
)
where npv stands for 'net pay value' and d i represents the value of the ith feature out of total n features. • The most probable net pay value (i.e. npv with the highest P(npv) value) is selected as the productivity index of corresponding depth, determined by the Bayesian classifier.

Fuzzy fusion

Sugeno integral is a famous nonlinear operator to fuse multisensory data. Fuzzy integrals have shown great results in fusing classifiers (Kuncheva 2004). The algorithm used here to apply Sugeno integral is presented below.

• Sort the components of the input vector (x) from the highest to the lowest value i.e.

x = [x 1 , x 2 , . . . , x n ] → x s = [x s 1 , x s 2 , . . . , x s n ] that : x s 1 > x s 2 > . . . > x s n .
• Arrange fuzzy densities of the components according to the sorted input vector (x s ), i.e.

g s = [g s 1 , g s 2 , . . . , g s n ]. • Calculateλ > -1 using relation (3): λ + 1 = n t=1
(1 + λg s t ).

(3)

• Set g(1) = g s 1 , and for t = 2 to n, g(i) should be calculated recursively using equation 4):

g(t ) = g s t + g(t -1) + λg s t g(t -1). ( 4
)
• Calculate the output of Sugeno integral using equation ( 5):

μ(x) = max t=1:n {min{x s t , g(t )}}.
(5)

Performance assessment

Confusion matrix and classification correctness rate (CCR) are used to evaluate the correctness of each classifier. The confusion matrix shows how input vectors are classified. The element (i, j) is the number of vectors that are classified in the ith class while belonging to the jth real class (Theodoridis and Koutroumbas 2003). The CCR is the ratio of the summation of corrected classified elements (usually on the trace of the confusion matrix) to the number of classes.

Results

Output of the conventional procedure

By using the conventional method, productive zones are determined crisply (figure 3). Cut-offs used to classify pay zones from non-pays are extracted from NIOC reports. For well no 1, 3, 4, 5 and 6, horizons with porosity values higher than 5%, water saturation of less than 40% and shale volume percentage of less than 25% are considered as pay zones. However for well no 2, the cut-offs of porosity and water saturation are considered as 6% and 50%, respectively, while the cut-off of shale volume is the same as in other wells.

As shown in figure 3, the whole Sarvak interval is divided into five different zones. This classification is the result of visual assessment. The upper half of Sarvak is more productive than the lower half, and zone number 2 is the most productive zone. Well no 4 and 5 are not incorporated in this figure because of inconsistency of their classification with the other wells. This inconsistency is probably due to imprecise classification of the Sarvak interval by the conventional procedure because the CCR value of well no 4 is considerably low in comparison to the other wells (figure 7(a)).

Output of the diffusivity equation-based methodology

For applying this methodology, in the first place, porosity and permeability values are estimated by an MLP-structured ANN, and water saturation is calculated by the Archie equation [START_REF] Johnson | Well Logging in Nontechnical Language[END_REF]. As it is described below equation (1), Q10 the distance from the well axis (r) and time (t) are considered as 5 ft and 5000 h respectively to keep the error of approximation less than 0.25%. Then, the index of productivity is calculated using equation (1), and presented in figures 6 and 4.

Q 0 p(r,
It was impossible to calculate the confusion matrix or CCR values for the outputs of this methodology because these performance assessment indices are compatible with discrete outputs, while the output of this methodology is continuous.

To check the validity of the results, the outputs of well tests are presented in the figure next to the outputs of the methodology. It is observed that productivity identification is at least acceptable within 14 well test intervals out of total 16 well tests. It was impossible to apply this methodology on well no 5 and 6 due to lack of core report, and consequently lack of porosity and permeability data. Generally, the upper part of Sarvak is more productive than the lower part, and zone number 2 has the highest productivity index.

Outputs of Bayesian classifiers

For training Bayesian classifiers using the mentioned algorithm, all the classes should exist in the training datasets. Therefore, datasets which can be used for training are well no 3 and 4 (see table 1). Furthermore in this research, Bayesian is trained twice by these two wells; one trained classifier is named Bayesian-3 (trained in well no 3), the other is called Bayesian-4 (trained in well no 4). The input well logs for training are LLD and LLD/LLS. Then, the trained classifiers are applied on all the wells to determine productive zones (figure 5).

As is obvious, the outputs of Bayesian-3 and 4 differ, to some extent, from each other (figures 5(a), (b)). Bayesian-3 identifies productive zones more pessimistically, whereas the output of Bayesian-4 is more optimistic. In this figure, like the outputs of previous methodologies, the upper half of Sarvak is more productive than the lower part. In the output of Bayesian-3 (figure 5(a)), zone number 2 is the most productive one, and in the output of Bayesian-4 (figure 5(b)), zone number 1 (equivalent of zone number 2 in Bayesian-3) is the most productive among the other zones. Bayesian-4 is unable to predict the productive zones through well no 6 precisely. It can be inferred from the low CCR value (figure 7), and the inconsistency of its output with other wells (figure 5(b)). The most significant difference between the outputs of Bayesian-3 and 4 is a weak productive zone identified by Bayeisan-4 (zone 4 in figures 5 and 6(b)), while Bayesian-3 has not determined this zone. In the next part, the outputs of Bayesian-3 and 4 are fused, by a fuzzy operator, Sugeno integral, to reach a more reliable result.

Output of fuzzy fusion

The outputs of Bayesian-3 and 4 are integrated by Sugeno operator to generate a new zoning. This zoning is more similar to the output of Bayesian-4 rather than that of Bayesian-3.

In this case, it seems that Sugeno integral has fused Bayesian classifiers optimistically. Incorrect classification of well no 6 is the worst disadvantage of this fusion-based zoning. However, higher certainty of this zoning in comparison to Bayesianbased zonings is the biggest advantage of this fuzzy-based zoning. The zoning is carried out by the means of visual evaluation, which is illustrated in figure 6.

Final integrated model for the productivity of the Sarvak formation

Based on the introduced methodologies, five zoning models for the productivity of the Sarvak reservoir are generated (figures 3-6). The differences, which exist between these figures, raise a considerable uncertainty in identifying productive zones of this reservoir rock. To solve this problem, these zonings are integrated by the technique of majority voting. In this integration, CCR values, calculated for each methodology, are used (figure 7). Sugeno integral are precise. However, due to CCRf values (figure 7(b)), the average of precision in Sugeno integral is not only the highest but also more constant well to well.

In the final zoning (integrating the five mentioned classifications), Sarvak formation is generally divided into two major parts: productive interval, which is approximately shallower than 3000 m (zones 1 to 3), and non-productive interval, which is almost deeper than 3000 m (zones 4 to 7). In a finer identification, this reservoir rock is divided into seven different classes (figure 8) from zone 1 to 7, in a descending order. Zone number 2 has the highest production rate, which makes it the most suitable horizon for perforation. Zones 1 and 7 are in the second priority and zone numbers 3 and 5 are in the third priority for production. Finally, zones 4 and 6 are not classified as oil producing zones. These horizons produce water instead.

Discussion

In this part, the final zoning is compared to previous sedimentological facts to check the validity of the proposed zoning. Sarvak is in three different depositional environments. Upper Sarvak belongs to fore reef with good reservoir quality; middle part belongs to main reef, which has the best reservoir quality, consequently very good productivity; and bottom Sarvak is a back reefal facies, related to low-energy environment with bad reservoir quality (Bashari 2007a). This sequence is comparable to the final results of this work too. Zone 1 with good reservoir quality is probably representative of fore reef deposition; zone 2 with the highest production rate is most likely related to the main reef with the best reservoir quality, and the productivity of other zones is similar to back reef i.e. bad reservoir quality. So, the proposed model is compatible with geological facts.

There are two zones which do not have lateral continuity in the southern part of the field (zones 3 and 7 in figure 8). It should be noted that these zones are not sedimentary zones; in fact, they are productive zones, which represent the production rate of that interval. This lateral variation may be due to different depositional environment or diagenetic processes. By the way, the effect of diagenetic processes is more probable because of two reasons: the first reason is that the variation in depositional environment rarely occurs within the field scale; on the other hand Sarvak has been under various diagenetic events such as dissolution, dolomitization, etc, affecting the porosity and permeability of the reservoir, consequently the rate of liquid production of the rock (Bashari 2007c, Taghavi et al 2007, Zamani et al 2010).

Conclusions

The conclusions of this paper can be categorized into two parts: methodological and geological. Methodological conclusions are the advantages and disadvantages of the introduced methodologies in identifying productive zones. Geological conclusions are those scientific facts and interpretations, inferred about the Sarvak reservoir in the field under study. The compatibility of the final method with geological facts, presented under discussion is the most valuable advantage of this work that verifies the application of the proposed methods in determining productive zones of Sarvak reservoir in this field.

Methodological conclusions. All the four newly developed methodologies identify productive zones more precisely than the conventional procedure. Besides, these proposed methods are capable of determining the productive zones fuzzily, while the conventional procedure classifies pay zones crisply. Due to CCRf values, Sugeno integral is the most precise fuzzy classifier among all the methodologies. In addition, the generalization ability of Sugeno operator is higher than that of other methodologies due to more constant CCRf values in different wells. At last, the continuous output of the diffusivity equation-based method is its obvious advantage among others.

Geological conclusions. On a large scale, Sarvak formation is divided into two major parts: productive interval, which is shallower than about 3000 m, and non-productive interval, which is approximately deeper than 3000 m. However on a finer scale, this reservoir rock is classified into seven zones and that zone number 2 has the highest production rate. In addition, the lateral variations of productive zones (3 and 7) are probably due to diagenetic processes, which happened on this giant reservoir. Page 9 Q11 Author: Please check the details for any journal references that do not have a blue link as they may contain some incorrect information. Pale purple links are used for references to arXiv e-prints.

Q13

Q12

Author: References 'Bashari (2007a)' and'Bashari (2007b)' appear to be the same. Please check.

Q13

Author: Please check whether the publisher name is okay as set in reference [START_REF] Johnson | Well Logging in Nontechnical Language[END_REF]

Figure 1 .

 1 Figure 1. Locationof Abadan Plain is illustrated in Zagros. Modified after Sherkati and Letouzey (2004) and Rajabi et al (2010).

Figure 4 .

 4 Figure 4. Output of the diffusivity-based methodology (index of productivity: ( Q 0 p(r,t ) )) in determining the productive zones of Sarvak formation.

Figure 5 .

 5 Figure 5. Outputs of Bayesian classifiers in determining the productive zones of Sarvak formation. White colour means npv = 3, orange means npv = 2 and black is npv = 1.

Figure 6 .Figure 7 .

 67 Figure 6. Output of Sugeno integral (fusing Bayesian-3 and 4) in determining the productive zones of Sarvak formation. White colour means npv = 3, orange means npv = 2 and black is npv = 1.

Figure 8 .

 8 Figure 8. Final zoning for the productivity of Sarvak formation, resulted from integrating the outputs of different classification methods.

Table 1 .

 1 Summary of datasets, available for this work. 'npv' stands for 'net pay value'. npv = 1 shows that the well test interval is not oil producing, npv = 2 reveals that the tested interval is oil producing but less than 1500 barrels per day and npv = 3 means that the well test result shows oil production more than 1500 barrels per day.

	Well 1 Well 2 Well 3 Well 4 Well 5 Well 6

•

  Intervals of well tests are marked by the production rate: if the interval does not produce oil, net pay value (npv) is defined as 1, if the corresponding depth produces oil less than 1500 barrels per day, then the net pay value equals 2 and finally if that interval produces oil more than 1500 barrels per day, then the net pay value is assumed as 3 in the corresponding well test interval.• Data of each mentioned class are divided into two parts: 70% for training Bayesian and 30% for testing the precision of trained Bayesian. • Probability distribution function of each feature in any of the three classes is calculated by training data.

	• For each test data (depth), the probability of belonging to
	each possible net pay value (npv = 1, 2 or 3) is calculated
	using equation (2):

  t ) Output of the conventional procedure in distinguishing net pay zones of Sarvak formation.

									W6
									2634
	W3	2697	W1	2726		W2
	2768
	3000	3000		3000				3000
									3270
		3334
	3382			3391				non productive N productive
	Figure 3. W3 w e ll t e s t id e n t if ie d p r o d u c t iv it y		2726 Perfect Acceptable Acceptable	W1			Acceptable 2726	W2	2626	W4
	2768		Acceptable					Perfect
	Perfect
			Acceptable
			Acceptable
	Acceptable
								Perfect
	3000 Acceptable		3000						3000	3000
	Acceptable
	Perfect		Perfect
									3271
			3334
	3382	N	highest production rate no production rate	0 2 4	x 10	-3	productivity index	3391

  The role of cut-offs in integrated reservoir studies SPE Reservoir Eval. Eng. 8 276-90 Zamani Z, Mirshahani M, Ghofrani E, Bayranvand B and Namati M 2010 Changes in reservoir quality determined from the diagenetic evolution of cretaceous carbonates: case study from Sarvak and Ilam Formations in the Lurestan Area, SW Iran 14th Int. Oil, Gas and Petrochemical Congress (Tehran) Please be aware that the colour figures in this article will only appear in colour in the Web version. If you require colour in the printed journal and have not previously arranged it, please contact the Production Editor now. Q3 Author: 'Khatiyeh' or 'Khatiyah'? Please check throughout.Q4Author: Please check whether the edits made in the sentence 'Bashari has completed another investigation . . . ' retain your intended sense. Please check whether the edits made in the sentence 'However it is claimed by Bashari that . . . ' retain your intended sense.Q7Author: Please check whether the edits made in the sentence 'The most important result of this research is . . . ' retain your intended sense. Please check whether the change 'national' from 'nationality' as given elsewhere in the text is correct here. Please check whether the edits made in the sentence 'The algorithm used here is a parametric Bayesian . . . ' retain your intended sense.Q10Author: The sense of the text 'As it is described below equation (1)' is unclear. Please explain what is meant.

	Kuncheva L L 2004 Combining Pattern Classifiers, Methods and Algorithms (Hoboken, NJ: Wiley-Interscience) Masoudi P, Tokhmechi B, Zahedi A and Jafari M A 2011 Developing a method for identification of net zones using log data and diffusivity equation J. Min. Environ. 2 accepted in August 2011 Piryaei A, Reijmer J J G, Borgomano J and van Buchem F S P 2011 Late Cretaceous tectonic and sedimentary evolution of the Bandar Abbas area, Fars Region, Southern Iran J. Pet. Geol. 34 157-80 Piryaei A, Reijmer J J G, van Buchem F S P, Yazdi-Moghadam M, Sadouni J and Danelian T 2010 The influence of Late Cretaceous tectonic processes on sedimentation patterns along the northeastern Arabian plate margin (Fars Province, SW Iran) Geological Society, London, Special Publications vol 330 pp 211-51 Rajabi M, Sherkati S, Bohloli B and Tingay M 2010 Subsurface fracture analysis and determination of in-situ stress direction using FMI logs: an example from the Santonian carbonates (Ilam Formation) in the Abadan Plain, Iran Tectonophysics 492 192-200 Shahvar M B, Kharrat R and Mahdavi R 2009 Incorporating fuzzy logic and artificial neural networks for building a hydraulic unit-based model for permeability prediction of a heterogeneous carbonate reservoir Int. Petroleum Technology Conf. (Doha, Qatar) Sherkati S and Letouzey J 2004 Variation of structural style and basin evolution in the central Zagros (Izeh zone and Dezful Embayment), Iran Mar. Pet. Geol. 21 535-54 Taghavi A A, Mørk A and Kazemzadeh E 2007 Flow unit classification for geological modelling of a heterogeneous carbonate reservoir: cretaceous Sarvak formation, Dehluran Field, SW Iran J. Pet. Geol. 30 129-46 Theodoridis S and Koutroumbas K 2003 Pattern Recognition (Elsevier) pp 489-544 Worthington P F 2010 Net pay-what is it? What does it do? How do we quantify it? How do we use it? SPE Reservoir Eval. Eng. 13 812-22 Worthington P F and Cosentino L 2005 QUERIES Page 1 Q1 Author: Keyword 'data fusion' is not cited in the text. Please check. Q2 Author: Page 2 Q5 Author: Please check whether the captions to figures 1 and 2 are okay as set. Q6 Author: Page 3 Q8 Author: Page 4	Q14 Q15 Q16 Q17
	Q9
	Author:

  Please provide page range in references 'Shahvar et al (2009)' and 'Zamani et al (2010)'.Q17Author: Please provide the place of publication in reference 'Theodoridis, Koutroumbas (2003)'. Also check whether the publisher name is okay as set.

	'.
	Q14
	Author: Please update reference 'Masoudi et al (2011)'.
	Q15
	Author: Please check whether reference 'Piryaei et al (2010)'
	is okay as set.
	Q16
	Author:

Acknowledgments

The authors wish to acknowledge Exploration Directorate of National Iranian Oil Company (NIOC) for permission to access company data. The authors thank Faraz Farbakhsh and Rouhallah Keshvardust for their helpful literal review on this paper. They also thank Reze Farajzadeh for his geological advices.


Pedram Masoudi
email: masoudip@ut.ac.ir

Yousef Asgarinezhad

Behzad Tokhmechi

Feature selection for reservoir characterisation by Bayesian Network

Keywords: feature conditioning, porosity, permeability, fracture, vug, net pay

The more accurate feature identification, the more precise reservoir characterisation.

Porosity, permeability and other rock properties could be estimated and classified by analytical and intelligent methods. Feature selection, plays a vital role in the process of identification. In this work, two goals are followed: first, developing Bayesian Network, K2 algorithm, as a complementary means (not an alternative) to find interrelationships of petrophysical parameters.

Second, feature conditioning for estimating porosity and permeability, vug and fracture detection, and net pay determination. Due to the results, bulk density log is introduced as the most important feature for characterising the reservoir because it is found useful for identifying all the studied reservoir features.

Introduction

Manuscript Click here to download Manuscript: BN,ver8.doc The concept of Bayesian Network (BN) was firstly developed in the fields of electrical and computer engineering. [START_REF] Pearl | Fusion, propagation, and structuring in belief networks[END_REF] and [START_REF] Cooper | A Bayesian method for the induction of probabilistic networks from data[END_REF] are of pioneers in Bayesian Network (BN) who defined this concept, and introduced the methodology clearly and applicably at the time. Later on, this methodology was used in a wide range of science and technology. [START_REF] Doguc | A generic method for estimating system reliability using Bayesian networks[END_REF] utilized BN in estimating system reliability.

Khor et al constructed three different types of BN classifiers in detecting network attacks; and by comparing the results, they concluded that these three types are well equivalent in performance [START_REF] Khor | From feature selection to building of Bayesian classifiers: A network intrusion detection perspective[END_REF]. BN is also used in some other fields like forecasting price in stock market [START_REF] Zuo | Stock price forecast using Bayesian network[END_REF]. It is some years that BN has been entered in geoscience studies. Based on the records of Scopus database; among all fields of earth science, remote sensing benefits from BN the most.

In petroleum industry, BN is used to assess situations and conditions probabilistically, e.g. in downstream it is used in circulation monitoring system [START_REF] Mansure | A Probabilistic Reasoning Tool for Circulation Monitoring Based on Flow Measurements[END_REF]; safety instrumentation and risk reduction at wellsite [START_REF] Kannan | Bayesian Networks: Application in safety instrumentation and risk reduction[END_REF]; identifying candidate wells for gel-polymer treatment [START_REF] Ghoraishy | Application of Bayesian networks for predicting the performance of gel-treated wells in the arbuckle formation, Kansas. 16th SPE/DOE Improved Oil Recovery Symposium 2008[END_REF]; drilling industry [START_REF] Al-Yami | Underbalanced Drilling Expert System Development[END_REF][START_REF] Al-Yami | Drilling Expert System for the Optimal Design and Execution of Successful Cementing Practices[END_REF][START_REF] Rajaieyamchee | Real time decision support in drilling operations using Bayesian Decision Networks[END_REF]; production issues and history matching [START_REF] Abdollahzadeh | Estimation of Distribution Algorithms Applied to History Matching[END_REF][START_REF] Hermann | Water Production Surveillance Workflow using Neural Network and Bayesian Network Technology: A Case Study of Bongkot North Field, Thailand[END_REF][START_REF] Khaz'ali | Bayesian network -A new probabilistic method for petroleum reservoir production prediction and history matching[END_REF]; completion [START_REF] Al-Yami | Expert System for the Optimal Design and Execution of Successful Completion Practices Using Artificial Bayesian Intelligence[END_REF]; and Enhanced Oil Recovery (EOR) [START_REF] Zerafat | Bayesian Network Analysis as a Tool for Efficient EOR Screening[END_REF] There are some publications of application of BN in upstream, specifically in basin analysis from economical evaluation of prospects [START_REF] Van Wees | A Bayesian belief network approach for assessing the impact of exploration prospect interdependency: An application to predict gas discoveries in the Netherlands[END_REF] to studying dependency relationships between geological features [START_REF] Martinelli | Bayesian networks for prospect analysis in the North Sea[END_REF][START_REF] Martinelli | Building Bayesian networks from basin-modelling scenarios for improved geological decision making[END_REF][START_REF] Rasheva | A new and improved approach for geological dependency evaluation for multiple-prospect exploration[END_REF]. In addition, there are two recently published papers in the upstream that have used BN in identifying effective logs, i.e. feature selection for determining productive zones through oil wells. Due to the results of one of articles, the ratio of LLD to LLS and individually LLD are the most effective raw features for detecting productive zones through oil wells (Masoudi et al., 2012c). Based on the results of the latter, porosity and water saturation are the most important extracted features for evaluating productive zones (Masoudi et al., 2012a).

It is worthy to mention that feature selection/ extraction is a basic and important stage in the process of identification [START_REF] Russo | Fuzzy methods for multisensor data fusion. Instrumentation and Measurement[END_REF]. It is not a good idea to consider all available information as input parameters. In another words, redundant information or duplications should be detected and removed from dataset [START_REF] Bleiholder | Data Fusion[END_REF].

The mentioned literature review reveals that newly developed concept of BN in petroleum industry is gradually going to become more and more applicable and renowned in exploratory investigations. As BN is a powerful tool to identify causal relationships between different features and phenomena, we have utilized it as a means to select effective petrophysical features for reservoir identification. The proposed procedure is based on correlation and dependency relations between reservoir properties and petrophysical parameters. In fact, we think that the deeper and the more precise understanding of interrelations and causations between parameters, the more effective feature selection, which plays an important role in success of any identification procedure; i.e. estimation, classification or clustering. Therefore, in this paper, authors follow two aims; the first one is developing the concept of dependency and Bayesian Network as an intelligent methodology for finding causality relationships and feature selection in petrophysical assessments, which is the novelty of this article. Second goal is introducing useful petrophysical parameters for identifying some reservoir properties (porosity, permeability, open fractures, vuggy porosity and net pay) within oil wells, which is a practical aid for petrophysiscists and geoscientists in their studies.

To do so, a brief review on a famous feature selection criterion, correlation coefficient, is presented following introducing available datasets; then, concept of "dependency" and methodology of "Bayesian Network" are added to make respected readers familiar with the concept and methodology. Thereafter, generated BNs and their outputs in various aspects of reservoir characterisation (estimating porosity and permeability, vug and fracture detection, and net pay determination) are included, followed by discussion and conclusion.

Datasets

In this work, petrophysical datasets of three Iranian oil-fields in Zagros Region have been studied. For the sake of confidentiality of data and information in National Iranian Oil Company (NIOC), the names of oil-fields under study (F1, F2 and F3) have not been enclosed but their approximate locations are indicated on Fig. 1. (Masoudi et al., 2012b;[START_REF] Rajabi | Subsurface fracture analysis and determination of in-situ stress direction using FMI logs: An example from the Santonian carbonates (Ilam Formation) in the Abadan Plain, Iran[END_REF][START_REF] Sherkati | Variation of structural style and basin evolution in the central Zagros (Izeh zone and Dezful Embayment), Iran[END_REF].

F1 is a giant field in Abadan Plain with North-South trend that has been used for evaluating net pay zones and estimating porosity and permeability. In this field, Sarvak Formation (Albian to Turonian) in six exploratory wells is studied. For fracture detection, one oil well in another giant oil-field (F2) is chosen. F2 is a Northwest-Southeast anticline in northern side of Kazerun Fault in South Dezful Area, very close to Izeh Zone. The reason why this field is selected for fracture study is availability of interpreted image logs and fullest of petrophysical data. For vug detection, a relatively small-sized anticline-shaped field (F3) in central Lurestan Area is selected. Access of authors to studied core reports is the reason for selecting this field to find causal relationship of vuggy porosity with petrophysical data. F3 has the same trend as F2, and like F1, investigation is fulfilled again within Sarvak Formation. Whereas the approximate locations could be seen in Fig. 1; summary of data and the purpose of choosing these three fields is summarized in Table 1.

Summary of available data and information in F1 are shown in Table 2. CGR, DT, NPHI, RHOB, LLD, LLS and MSFL are common well logs in all six wells; therefore, in order to incorporate maximum number of wells, other well logs are not included in this study. In addition, because there is no core data in well 6, this well is exempted from porosity-permeability study.

Also, due to lack of well test data in well 5, this well is exempted from net pay investigation.

In each of F2 and F3 fields, only one well is available. Available data in F2 are CGR, NPHI, DT, PEF, RHOB and SGR well logs, and interpretation of open fractures on image log; whereas available data in F3 are GR, Cali, RHOB, DRHO, NPHI, DT and LLD well logs, and observed vuggy porosity in cores.

A Simple Review of the Correlation Coefficient

Correlation coefficient is a well-known factor, measuring correlation (mutual relationship) between two different variables. There are different standpoints for calculating correlation coefficient: algebraic, geometric, and trigonometric. Pearson product-moment correlation coefficient is the most well-known formula for calculating correlation coefficient of two variables from algebraic viewpoint [START_REF] Lee Rodgers | Thirteen ways to look at the correlation coefficient[END_REF]. Fig. 2 shows two correlated variables; i.e. b changes when a changes in the same or reverse direction (variables in Fig. 2 are correlated in the same direction). Although correlation coefficient is a very valuable and important factor for understanding interrelations of variables, there are some insufficiencies in using it [START_REF] Bobko | A Review of the Correlation Coefficient and Its Properties, Correlation and Regression[END_REF]. does not necessarily mean that the number of children has risen too (Whereas Population Growth Rate is positive in developing countries, it is very close to zero or even negative in developed countries, and is not directly related to number of fathers or adults). But when the number of kids rises, you are 100% sure that the number of fathers (adults) has risen; because every kid needs a father to be born but fathers do not need their children for existing! Therefore, there is no mutual relation or correlation between number of fathers and children; however one of them is dependent on the other (directional relation).

Another example for insufficiency of correlation coefficient in showing interrelation of two variables, revealed in Fig. 3. If 100 people are asked to climb Damavand Mountain (highest peak in Iran with elevation of 5610 meters above the geoid), and plot the height of which they have reached against their ages, the plot would be like in Fig. 3. In fact, the acquired data is distributed between two envelopes that show possible and certain accessible heights for each age.

The shapes of these two envelopes are similar to an inverse "V", because teenagers and olds are not able in reaching high heights, whereas young people and middle-aged can even reach the peak. Correlation coefficient of dataset of this plot is something close to zero but without any shadow of doubt, there is a relationship between maximum accessible height and age in mountain climbing, while they are not correlated due to correlation coefficient.

Fig. 3. Dependency between age and accessible height for a man.

A practical example of this inefficiency in petrophysics could be found on cross-plot of Calliper-Vug and RHOB-Vug on Fig. 8. In these two cross-plots, there is no mutual relationship between two plotted variables as in Fig. 3, later we show that both Calliper and RHOB are important features for vug detection.

What Is Dependency?

This work introduces dependency between variables as an alternative for finding related variables, especially when there is no mutual relationship like relationship of father and child. Now, what is dependency? Each field has its own definition of dependency, and they are close to each other. Oxford dictionary states that "dependence" means "the state of relying on or being controlled by someone or something else" (OxfordDictionaries, 2010). From mathematical viewpoint, probabilistic is a means for evaluating dependency of variables on each other.

Probabilistically, two variables are called independence, when the joint probability of them is equal to product of their own probability [START_REF] Olofsson | Probability, statistics, and stochastic processes[END_REF]:

(1)

In the above equation, A and B are independent sets. Bayesian Probability is theory of studying conditional probability between two or more dependent variables, which uses Bayes rule for calculating evidential probabilities [START_REF] Duda | Pattern Classification[END_REF]. Bayesian Network, which is introduced here in order to find out dependency relation between petrophysical variables, is mainly based on Bayesian theory of conditional probabilities.

Methodology

Bayesian Network

Consider five measured variables (named a 1 to a 6 ) that are supposed to be effective on another unknown variable, called b, and each of these seven parameters can admit four different states. Therefore there are four powered seven (i.e. 16'384) possible states, and it is not only hard to compute and consider all these states (an NP-hard problem), but also impossible in the case of lack of complete and comprehensive dataset of records.

Bayesian Network is a directed acyclic graph that nodes represent variables and edges show direct dependencies between the linked variables. Now, suppose that dependency relationships between those previously mentioned seven variables can be represented as in Fig. 4. Based on this graph, the variable b is only dependent on the variables a 3 , a 4 and a 5 . Also, it is simple to formulate probability of b as:

(2)

(3) That P(x) is probability of occurrence of x, is joint probability, i.e. probability of occurrence of x and y simultaneously, is probability of x, considering y, i.e. conditional probability. The first equation is inferred from independency of the variable a 5 and set of variables a 3 and a 4 . For better understanding of equation 3, respected readers are referred to [START_REF] Pearl | Fusion, propagation, and structuring in belief networks[END_REF]. There are two methodologies for constructing BNs: constraint-based methods and scorebased methods [START_REF] Lauría | An Information-Geometric Approach to Learning Bayesian Network Topologies from Data[END_REF]. The former is used in cases that user is confident about the causal relationships between variables. For instance Total Organic Carbon (TOC) is a prerequisite for oil generation, and none of specialists believe that oil could be generated without having some least amount of TOC [START_REF] Al-Ameri | Petroleum system analysis of the Mishrif reservoir in the Ratawi, Zubair, North and South Rumaila oil fields, southern Iraq[END_REF]. Furthermore there is a dependency relation between TOC and oil generation. In fact, constraint-based methods are judgmental methods that an expert is responsible for [START_REF] Martinelli | Bayesian networks for prospect analysis in the North Sea[END_REF].

In some cases, it is difficult or even impossible for a user to determine dependency relations between variables. In these cases, data-driven approaches are used to find the most probable state of dependency between each pair of variables. In score-based methods, a calculated score is set as a criterion to find the dependency relation between two variables. For using score-based methods, two elements should be specified: search procedure and scoring metric. Scoring function should be associated with probability of a candidate directed acyclic graph, and search procedure is considered as an optimization problem. Greedy hill-climbing algorithm, K2 algorithm, simulated annealing optimization, Monte Carlo are some of those score-based methods, known as heuristic approaches to construct a BN [START_REF] Cooper | A Bayesian method for the induction of probabilistic networks from data[END_REF][START_REF] Lauría | An Information-Geometric Approach to Learning Bayesian Network Topologies from Data[END_REF][START_REF] Niedermayer | An introduction to Bayesian networks and their contemporary applications[END_REF]. In this work, K2 algorithm is used in order to construct BNs.

K2 Algorithm

K2 algorithm is a score-based method for constructing BNs, although it is not completely free of constraint. Two constraints should be considered prior to running K2 algorithm. The first one is to state the maximum possible parents that each node can have; the other is providing a true initial order of variables by user that variables are not dependent in reverse order. E.g. if the initial order of (a 1 , a 2 , a 3 , … , a n ) is provided by user, a j could be dependent on a i (i<j); though a i cannot be dependent on a j but if a i and a j are mutually correlated. After considering these two constraints, following algorithm should be applied on the dataset [START_REF] Doguc | A generic method for estimating system reliability using Bayesian networks[END_REF]:

Algorithm K2(T,u):

Input: dataset of observations, T, and maximum possible number of parents for each node, u.

Output: Bayesian Network, BN.

(1) For each variable in input dataset, T

(1-1) Create node Ai as i-th variable, and add it to BN

(1-2) Create an empty set as set of parents (Pa i ) of A i

(1-3) Calculate by:

(4)

Where, is number of times that A i and A j are in a specific state of k. d i is number of states of A i . Finally, q i is number of possible parents, i.e. .

(1-4) While number of elements of Pa i is not larger than u:

(1-4-1) Assume X z as a parental node of X i

(1-4-2) Calculate

(1-4-3) If the score of is larger than , fix A z as a permanent node of A i , otherwise, remove it from parental set (Pa i ).

(2) Return BN

Results

In this section, correlation coefficient and Bayesian Network are utilized to find out effective features for reservoir identification. In the first part, important features for porosity and permeability estimation are determined, and in the second, third and fourth parts, useful features for fracture and vug detection and net pay assessment are determined respectively.

Causal Relationships in Porosity and Permeability Estimation

Estimating porosity and permeability of reservoir rocks is very essential in reservoir characterization, static and dynamic modelling. There are many investigations about estimating porosity and permeability of reservoirs. The input features, used for estimating these two parameters have not remained unchanged during time. Table 3 shows different datasets, introduced for porosity and permeability estimation in chronological order.

In order to find out, appropriate features for estimating porosity and permeability, scatter plots of all logs, and porosity and permeability of core analysis are plotted on Fig. 5 Due to criterion of correlation coefficient, porosity is relatively well-correlated to Log(Perm) and RHOB. Based on constructed BN, porosity is related to Log(Perm) (child) and RHOB (parent) in the first order. Therefore, both criteria, unanimously agree that the first two features, related to porosity are Log(Perm) and RHOB. Again due to correlation coefficient, Log(Perm) is relatively well-correlated with porosity and LLD. BN approves that porosity is the most influential factor for estimating Log(Perm), whereas defines RHOB in the second stage.

It is relatively easy to come on an agreement in feature selection for porosity because priorities of both criteria are very close to each other. Three features of RHOB, DT and NPHI are the most related features to porosity; then, LLD and CGR. For Log(Perm), it is a bit tricky. The two mostly related features for permeability estimation are porosity and RHOB. After these two, LLD, NPHI and DT could be named.

Causal Relationships for Fracture Detection

Studying fractures is much more complex than porosity and permeability due to wild nature of fractures and variety of fracture geometry. Image logs are main means to identify and characterize fractures [START_REF] Ja'fari | Fracture density estimation from petrophysical log data using the adaptive neuro-fuzzy inference system[END_REF], although in case of lack of this information source, traditional well logs are used (Table 4). Surprisingly, based on both criteria, priorities of features for fracture identification are the same, except in the places of CGR and PEF.

Causal Relationships for Vug Detection

Like fractures, image logs are most reliable tools for vug detection. It is worthy to mention that vug pores are visible in cores likewise, whereas fractures could not be studied in cores due to low core recovery within fractured intervals. In addition, vug pores are not investigated as much as fractures so far (Table 4), due to their relatively low importance, comparing to fractures. By the way, for selecting appropriate features for vug detection, scatter Calliper log is the most important feature for vug detection (Fig. 9 (a)); NPHI, DT, RHOB and GR are other important features in order.

Causal Relationships for Net Pay Determination

Determining productive zones is a very critical stage in static reservoir modelling.

Petrophysical net pay determination is usually done by cut-off method. Some of utilized features in literature are included in Table 5.

For net pay detection, production rate, derived from production test is utilized as criteria of productivity: In order to satisfy the second goal of the paper, which is providing results of feature selection in order to benefit petrophysiscists of, selected features for reservoir characterization are summarized in Table 6. 1 st stage features are those features that have high priority due to both criteria; 2 nd stage features are those effective features that do not have the same importance as 1 st stage features. For fracture detection, there is another column named 3 rd stage features that are not as important as 2 nd stage features. It is worthy to mention that bulk density (RHOB) is the most frequent feature in this table; therefore, the most important log for reservoir characterization.

Tables 3, 4 and 5 are included in the current work in order to validate obtained results, i.e. proposed input features for reservoir study, Table 6. Comparing the current (Table 6) and previous works (Tables 3, 4 and5) reveals that selected features (Table 6) are reasonable for porosity, permeability, fracture, vug and net pay studies; furthermore, combining BN and correlation coefficient is a successful way for feature extraction in reservoir characterisation.

Conclusion

Although correlation coefficient is a very useful and easy to use criterion to find and quantify mutual relationships between different variables, there are some pitfalls when using it.

In this work, Bayesian Network is introduced as a complementary means (not an alternative) to find out dependency relations; therefore, finding causal relationships and feature selection in

Fig. 1 .

 1 Fig. 1. Map of Zagros Region and main faults of the region. Modified after (Masoudi et al.,

Fig. 2 .

 2 Fig. 2. The more correlated variables (here a and b), the closer dots to the dashed line.

Fig. 4 .

 4 Fig. 4. A Bayesian Network, showing dependency relation of seven variables of a 1 to a 6 and b.

Fig. 5 .

 5 Fig. 5. Relationships between features in Sarvak Formation of F1 for porosity and permeability estimation.

  Like before, to find out appropriate features for fracture study, scatter plots of all logs, and identified open fracture on image logs are shown on Fig. 6. Then, correlation coefficients between each pair of variables are calculated and are shown on each plot. Due to correlation coefficients, DT and SGR are the most important features for fracture identification in F2; then, CGR and RHOB.

Fig. 6 .FractureFig. 7 .

 67 Fig. 6. Correlation chart and histograms in Sarvak Formation of F2.

  plots of all logs, and observed vug pores are shown on Fig. 8. Then, correlation coefficients between each pair of variables are calculated and shown on each plot. Due to correlation coefficients, NPHI, DT and RHOB are the most important features for vug detection in F3.

Fig. 8 .Fig. 9 .

 89 Fig. 8. Correlation chart in and Sarvak Formation of F3.

Fig. 10 .

 10 Fig. 10. Relationships between features in Sarvak Formation of F1 for net pay detection.

Figure1

  Figure1Click here to download high resolution image

Table 1 .

 1 The table shows what reservoir properties are studied in which field

	Figure4

Table 2 .

 2 Summary of dataset of F1 oil field, available for evaluating causality relationships for assessing porosity, permeability and net pay zones

	Well 1 Well 2 Well 3 Well 4 Well 5 Well 6

Table 3 .

 3 Petrophysical parameters for porosity estimation in various references

	Table3
		Used Parameters for Estimation	Source
		NPHI, Density, Sonic, Resistivity	(Helle et al. 2001)
	Porosity	DT, GR, ILD, ILS, NPHI, RHOB	(Jalali Lichaei and Nabi Bidhendi 2006)
		CGR, DT, LLD, LLS, MSFL, NPHI, RHOB	(Masoudi et al. 2011b)
		NPHI, Density, Sonic, Resistivity	(Helle et al. 2001)
		Depth, DT, GR, ILD, ILS, RHOB, Sw, Porosity	(Jalali Lichaei and Nabi Bidhendi 2006)
		NMR	(Fethi et al. 2010; Timothy et al. 2008)
		SGR, CGR, RHOB, TNPH (Thermal Neutron Porosity),
	Permeability	Rs (medium resistivity), Rt (deep Resistivity), Rxo (shallow resistivity), DT, VCLAY (clay volume) Saturation, Gamma, Neutron, RHOB, PEF, DT, Resistivity	(Shahvar et al. 2009) (Saemi et al. 2007)
		Gamma, DT, Nphi, RHOB, LLD/LLS	(Ibrahim Sami and Adel 2010)
		NPHI, RHOB, DT, LLD, SGR, CGR	(Mehri 2010)
		CGR, DT, LLD, LLS, MSFL, NPHI, RHOB	(Masoudi et al. 2011b)

Table 4 .

 4 Petrophysical parameters for evaluating secondary porosity

	Table4
		Used Parameters for Identification	Source
		Water Saturation, GR	(Tokhmechi et al., 2009)
	Fracture	Calliper, DT, RHOB, PEF	(Tokhmchi et al., 2010)
		DT, RHOB, NPHI, Resistivity	(Ja'Fari et al., 2012)
	Vug	NPHI, DT, GR, Calliper NPHI, DT, GR, RHOB	(Asgarinezhad et al., 2011) (Asgari-Nezhad et al., 2012)

Table 5 .

 5 Petrophysical parameters for net pay determination in various references

	Table5
	Used Parameters for detection	Source
	Shale Volume, Porosity, Water Saturation	(Jensen and Menke 2006; Mahbaz et al. 2011; Worthington 2010)
	Permeability, Porosity, Viscosity, Compressibility	(Masoudi et al. 2011b)
	Porosity, Water Saturation, Shale Volume	(Masoudi et al. 2012a)
	Ratio of LLD to LLS and LLD	(Masoudi et al. 2012c)

Table 6 .

 6 Result of feature selection for porosity, permeability, fracture detection, vug detection and net pay determination due to correlation coefficient and Bayesian Network

Acknowledgment

The authors wish to acknowledge Exploration Directorate of National Iranian Oil Company (NIOC), for permission to publish scientific results. Also, special thanks to Mr.

Meisam Akbari Zare, a professional mountain climber for sharing his nice photo of Damavand Mountain with us. Finally, many thanks to Dr. Abdullah M. Al-Amri and Dr. Jassim M. Thabit,


Pedram Masoudi
email: masoudip@ut.ac.ir

Hossein Mohammad

Net Pay Determination by Dempster Rule of Combination, Case Study on Iranian Offshore Oil Fields

Keywords: production rate, net pay, uncertainty, well test, Dempster-Shafer Theory

Net pay detection is a key stage in reservoir characterization for several purposes: reserve estimation, reservoir modeling and simulation, production planning, etc. Determining productive zones always is simultaneous with some amount of uncertainty due to lack of enough data, insufficiency of knowledge and wild-nature of petroleum reservoirs. It becomes even more challenging in carbonates, because of their highly heterogeneous environment. Conventionally, evaluating net pays is done by applying petrophysical cut-offs on well-logs, which results in crisp classification of pay or non-pay zones.

In addition, cut-off-based method is developed in sandstones, and does not provide suitable results in carbonates at all. Proposed methodology of this work, Dempster-Shafer Theory, is a generalization of Bayesian Theory of conditional probabilities. Net pays are studied in two oil reservoirs by this theory; one of them is carbonate reservoir of Mishrif, the other is sandy Burgan reservoir. For validation, results are compared to well tests and output of conventional cut-off method. The advantages of using Dempster-Shafer Theory, comparing to conventional cut-off-based method in studying net pays is: to have a continuous fuzzy output, based on geological facts, with high generalization ability and more compatibility with well test data.

Introduction

We drill to produce; therefore, finding productive zones (net pays) backs to the very first drilled wells in the story of oil industry. The history of net pay determination could be categorized into three parts. In the first part, specialists used to determine net pays while drilling. Second part starts just after invention of logging tools. In this part, finding productive zones is based on indirect methods; i.e. petrophysicists used to "infer" where the oil is accumulated with appropriate reservoir pressure. Third part is followed by development of computer and computing systems in net pay determination. In this part, there is not a significant hardware progress but software are highly improved, based on acquired data, mostly logs. The methodologies vary from simple statistical assessments to complex data-driven methodologies. The below chart is summary of utilized tools or methods for net pay determination: Direct Methods [START_REF] Millikan | Use of Gas Meters for Determination of Pay Strata in Oil Sands[END_REF] Gas-meter [START_REF] Connell | Development of a Fluorimetric Method for Detection of Pay Zones During Drilling With Invert Muds[END_REF] Flourimetric [START_REF] Chaudhry | Introduction, Oil Well Testing Handbook[END_REF][START_REF] Matthews | Pressure Buildup and Flow Tests in Wells[END_REF] Drill Stem Test (DST) and other well testing methods [START_REF] Cooke-Yarborqugh | Reservoir Analysis By Wireline Formation Tester: Pressures, Permeabilities, Gradients And Net Pay[END_REF] Cased-hole Wireline Formation Tester New tool for well testing Part II: Infera ble Meth ods [START_REF] Snyder | A Review of the Concepts and Methodology of Determining "Net Pay[END_REF] Gamma and resistivity logs [START_REF] Flower | Use of Sonic-Shear-Wave/Resistivity Overlay as a Quick-Look Method for Identifying Potential Pay Zones in the Ohio (Devonian) Shale[END_REF] Resistivity and shear sonic logs [START_REF] Cooke-Yarborqugh | Reservoir Analysis By Wireline Formation Tester: Pressures, Permeabilities, Gradients And Net Pay[END_REF] Pressure meter New instrument as a quick-look [START_REF] Deakin | The Integration of Petrophysical Data for the Evaluation of Low Contrast Pay[END_REF] Conventional logs, capillary pressure of cores and interpreted image logs Gas reservoirs [START_REF] Worthington | Recognition and evaluation of low-resistivity pay[END_REF] Gamma, caliper and self-potential logs recognizing low-resistivity pays, tight sands [START_REF] Nottale | Scale relativity and fractal space-time: applications to quantum physics, cosmology and chaotic systems[END_REF] Conventional logs and core data Application of rock typing in net pay detection Part III: Optimized and data-driven Methods [START_REF] Worthington | The Role of Cut-offs in Integrated Reservoir Studies[END_REF] Shale volume, porosity and water saturation Review paper on cut-off method [START_REF] Jensen | Some Statistical Issues in Selecting Porosity Cutoffs for Estimating Net Pay[END_REF] Cut-off estimation Statistical issue for estimating cut-off [START_REF] Worthington | The Application of Cutoffs in Integrated Reservoir Studies[END_REF] Dynamic cut-off ( Authors believe that the answer depends on the standpoint of which we are looking from. In general for petroleum industry, there are two viewpoints: "Engineering" and "Science". It is ideal for engineers to design the easiest and fastest pathway towards destination, meanwhile meeting required criteria, accuracy and precision. On the other hand, for scientists, it is essential to explore all possible ways, categorizing methodologies and achievements with less economic considerations, in order to benefit of findings, as is required by engineers. Hence, science is a support for engineering. In another word, science is a domain for engineering, and the more progressed the science is, the wider domain for engineers to act and build within. Therefore having multiple solutions for a single problem is power, and improves engineer's abilities, and confidence of their designing and decisions.

Specifically about the case of "Net Pay Detection", we should consider highly heterogeneous condition of reservoirs. Reservoir conditions vary a lot. This heterogeneity results in uncertainty of geo- In this situation, some methodologies might not be applicable in some conditions, e.g. some methodologies are designed for simpler sandy reservoirs, and not applicable in wild-nature carbonate reservoirs or vertically heterogeneous environment of laminated reservoirs. Hence, having multiple approaches makes industry capable to handle various conditions. In this work, two reservoir conditions are considered in order to check proposed DST methodology in two prominent reservoir environments: carbonate and clastic.

Datasets and Geologic Settings

Developed methodology is applied in two reservoirs: a carbonate and a sandy reservoir of two different oil fields. Both structures are of Iranian offshore oil fields in the Persian Gulf, and geological aspects of studied reservoirs are described as follow. Utilized data in both datasets consists of calculated porosity, shale content and water saturation, in addition to well test results for verification. Sandy Burgan reservoir, Middle Cretaceous (Albian era), is a Member of Kazhdumi Formation [START_REF] Bashari | Integrated 3D seismic and petrophysical data of the Sarvak formation in the Persian Gulf[END_REF][START_REF] Mahbaz | Optimization of reservoir cut-off parameters: a case study in SW Iran[END_REF]. As it is reported from the studied field (we render this field as "S" field due to studying a sandy reservoir through), Burgan Member is overlaid on unnamed clastics, after a sedimentation gap. Sequence stratigraphy of Burgan Member in an adjacent offshore oil field, Forouzan;

shows that Burgan consists of four sedimentation sequences that each one starts from shallow sandy facies, and ends with deeper shaly facies [START_REF] Honarmand | Lithofacies Analysis, Sedimentary Environment and Sequence Stratigraphy of Burgan Formation, North of the Persian Gulf[END_REF]. The available dataset in this oil field is comprised of five drilled wells.

Methodologies

Net pay determination is done by two methodologies: conventional and proposed. The conventional method is cut-off-based method, which provides a crisp output. The results of cut-off method are considered as a bench mark for assessing proposed method. The developed method is based on Dempster-Shafer Theory (DST). The methodologies are introduced as follow.

Conventional Cut-off-based Methodology

Conventional cut-off-based methods provides a crisp (0/1) output (Masoudi et al., 2012a;Masoudi et al., 2012b;Masoudi et al., 2012c;[START_REF] Masoudi | Developing a Method for Identification of Net Zones Using Log Data and Diffusivity Equation[END_REF][START_REF] Worthington | Net Pay-What Is It? What Does It Do? How Do We Quantify It? How Do We Use It[END_REF].Due to reports of Iranian Offshore Oil Company (IOOC), kaolinite is dominant shale mineral in Mishrif reservoir of "C" oil field; furthermore, kaolinite-porosity plot is produced in each well, in order to distinguish different reservoir facies. Then, borders of each reservoir facies is interpreted as cut-off value of porosity or kaolinite. Cut-off of water saturation is inferred from histograms of water saturation in each well, where two distributions are distinguishable. Figure 1 shows procedure of net pay determination by cut-off method.

Figure 1Procedure of net pay determination by cut-off methodology.

Dempster-Shafer Theory (DST)

"Dempster-Shafer Theory of Evidences" [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF], "Bayesian Theory of Conditional Probability" and "Possibility Theory" [START_REF] Negoita | Fuzzy sets as a basis for a theory of possibility[END_REF] are three major schools of thought in assessing Body of Evidences and reasoning. Dempster-Shafer Theory (DST) is generalization of Bayesian Theory, developed under the concept of "Subjective Probability" [START_REF] Shafer | Perspectives on the theory and practice of belief functions[END_REF]. Reputation of DST is in decision-making in uncertain situations [START_REF] Rakowsky | Fundamentals of the Dempster-Shafer theory and its applications to system safety and reliability modelling[END_REF]; like here in uncertain condition of petroleum reservoirs. In fact, DST is another mathematical approach for assessing uncertainty [START_REF] Sentz | Combination of Evidence in Dempster-Shafer Theory[END_REF].

Aminzadeh is the first researcher who has used DST in uncertainty assessment in petroleum exploration. It is reported in his research that by fusing opinions of two specialists about probability of oil accumulation in a reservoir, uncertainty is lowered significantly. In his sample, the first expert has claimed that the chance of oil accumulation in a prospect is at least 40% (i.e. 60% ignorance). The second expert has claimed that the probability of existence of water or oil is at least 60% (i.e. 40% ignorance). By using DST to combine these two opinions, it is inferred that chance of finding "oil" is between 40% and 64% (i.e. 64%-40%=24% ignorance); and chance of finding "water or oil" is from 36% to 60% (i.e. 60%-36%=24% ignorance). Therefore the ignorance range is lowered from 60% or 40% to 24%.

As in the above sample, output of DST is a range of probability, and boundaries are called "Plausibility" and "Belief", and could be calculated by "Plausibility Function" and "Belief Function"

respectively. There are some formula for calculating Belief and Plausibility but the easiest way is by using another function called "Mass Function", which should satisfy the below conditions:

𝑚 ∅ = 0 (3) 𝑚 𝐴 𝑖 2 𝑛 𝑖=1 = 1 (4)
Where, "m" is Mass Function; "A i " is a possible state out of totally "2 n " states; while "n" is number of elements/ variables in the system under investigation. Based on Mass Function, Belief and Plausibility Functions are defined as:

𝑏𝑒𝑙 𝐴 𝑖 = 𝑚 𝐵 𝑗 𝐵 𝑗 ⊆𝐴 𝑖 (5) 𝑝𝑙 𝐴 𝑖 = 𝑚 𝐵 𝑗 𝐴 𝑖 ∩𝐵 𝑗 ≠0 (6)
In fact, Belief Function is a pessimistic viewpoint about occurrence of "A i "; whereas, Plausibility Function is an optimistic standpoint. In mathematical language, the relationship between these functions and probability is:

𝑚 𝐴 𝑖 ≤ 𝑏𝑒𝑙 𝐴 𝑖 ≤ 𝑝 𝐴 𝑖 ≤ 𝑝𝑙 𝐴 𝑖 (7)
In the case of a multisensory system, Mass Function could be calculated for the whole system by using Dempster Rule of Combination. This rule (equations 8 and 9 are Dempster Rule of Combination for a two-sensory system) uses Mass Functions of sensors (m 1 and m 2 ) in order to produce a new Mass Function for decision-making about whole the system; in another word, for sensor fusion [START_REF] Liu | Classic Works of the Dempster-Shafer Theory of Belief Functions: An Introduction[END_REF][START_REF] Rakowsky | Fundamentals of the Dempster-Shafer theory and its applications to system safety and reliability modelling[END_REF]. It is essential to mention that Dempster Rule of Combination is a purely conjunctive (AND) operation; and mathematically, it is recommended to combine independent bodies of evidences, (m 1 and m 2 ) [START_REF] Sentz | Combination of Evidence in Dempster-Shafer Theory[END_REF]; however this constraint is not a big restriction for engineering use if output satisfies requirements.

𝑚 𝐴 = 𝑚 1 𝐵 𝑚 2 𝐶 |𝐵 ∩ 𝐶 = 𝐴 𝑚 1 𝐵 𝑚 2 𝐶 |𝐵 ∩ 𝐶 ≠ ∅ (8) 𝑚 𝐴 = 𝑚 1 𝐵 𝑚 2 𝐶 |𝐵 ∩ 𝐶 = 𝐴 1 -𝑚 1 𝐵 𝑚 2 𝐶 |𝐵 ∩ 𝐶 = ∅ (9)

Proposed DST-based Methodology

For proposed method, the same inputs of conventional cut-off method, i.e. porosity, shale content and water saturation are used to determine net pay zones. The authors selected the same input set forboth conventional and proposed methods just to examine functionality of the proposed method versus conventional method apart from the effect of input selection. Then, for each input parameter, the simplest Mass Function is created in a way that higher value shows good reservoir quality:

𝑚 𝑝𝑜𝑟𝑜𝑠𝑖𝑡𝑦 = 𝑝𝑜𝑟𝑜𝑠𝑖𝑡𝑦 (10) 𝑚 𝑠ℎ𝑎𝑙𝑒 = 1 -𝑠ℎ𝑎𝑙𝑒 (11)
𝑚 𝑤𝑎𝑡𝑒𝑟 𝑠𝑎𝑡𝑢𝑟𝑎𝑡𝑖𝑜𝑛 = 1 -𝑤𝑎𝑡𝑒𝑟 𝑠𝑎𝑡𝑢𝑟𝑎𝑡𝑖𝑜𝑛 (12)
Thereafter, Dempster Rule of Combination (equation 8 or 9) is used in order to fuse the predefined Mass Functions (10, 11 and 12), which we call it Mass Function of net pay (whole the system) or output of proposed DST-based method [START_REF] Masoudi | Identifying Net Pay Zones in Oil Wells by Artificial Neural Network and Dempster-Shafer Theories[END_REF]. Mathematically, Dempster Rule of Combination is invented for the case of integrating information of independent sources [START_REF] Sentz | Combination of Evidence in Dempster-Shafer Theory[END_REF]. In this work, fusing data sources are porosity, shale percentage and water saturation. In engineering applications, sometimes users do not take account this limitation for using the rule due to precise and accurate outputs. In another word, this constraint is not a rigid one for engineering applications. Here, in both reservoirs, water saturation does not show any correlation with other two variables; i.e. correlation coefficient is close to zero in all wells. In carbonate reservoir of Mishrif, there is a low rate of correlation. But in sandy Burgan reservoir, correlation coefficient is at maximum: around 0.6, which could be categorized as semi-dependent. Hence, in this study, it is logical to use Dempster Rule of Combination, mathematically.

Resulted Mass Function gives a value from zero to much lower than one for each horizon.

Summation of all Mass Function values in each well is equal to one due to definition of Mass Function.

Finally, summation of Mass Function over a fixed interval, gives a new function, called Belief function, which shows degree of certainty of that interval to be productive. By investigating effect of neighboring horizons on production, when producing from a perforated interval (vertical permeability and viscosity are required), Plausibility Function could be easily calculated by widening the interval of summation with respect to vertical permeability and viscosity values.

Evaluation Criterion

In order to evaluate the proposed methodology, the output is compared with conventional cutoff-based method. Then, closeness of outputs of both methods is compared with well-test results. From well-test data, oil rate is normalized (divided) by interval length of testing interval, which gives barrel oil production per day per meter ( 𝑏𝑏𝑙 𝑜𝑖𝑙 𝑑𝑎𝑦 .𝑚 ), and is plotted against output of both methodologies. In another words, oil rate is normalized by length of testing to remove the effect of testing length on oil rate.

Results and Discussion

Based on explained procedure in Figure 1, cut-offs of porosity, shale volume and water saturation are determined in each well of two available datasets. For carbonate reservoir of Mishrif, the range of porosity cut-off was from 10% to 20%. Please consider that, relatively whole the Mishrif is ahigh-quality reservoir (low-shale and high-porosity), and cut-offs of porosity and shale volume are only defined to distinguish the cleanest pay zones from whole the interval to have a pay-categorization through wells. Cut-off of shale volume was below 10%; and cut-off of water saturation was from 50% to 80% due to distribution variability of water saturation in each well.

In the sandy reservoir, defined cut-offs were closer to ordinary values because Burgan reservoir is not as clean as Mishrif. Porosity cut is defined from 4% to 11%; shale volume between 42% and 60%;

and finally water saturation from 35% to 61%.

Output of DST-based methodology, conventional method and testing results are plotted in ).

1. Both DST and cut-off-based methodologies are compatible in detecting non-pay zones. In another word, the results would be the same, if someone identifies non-productive zones by each procedure. Please note zero values in Figure 2 and3.

2. In geology, the borders are mostly transitional, not sharp. Although there is a sharp change in reservoir parameters in the contact of reservoir-cap rocks, changes of reservoir parameters within the reservoir interval is transitional. E.g. water saturation does not convert sharply from say 20% to 80%, it is gradual. Also the same is for porosity, shale volume and other reservoir parameters. Therefore, dividing the gross into pay and non-pay zones is unrealistic. Cut-off-based method gives a crisp (0-1) output, which is unrealistic within the reservoir rock, while DST-based method gives a more realistic, fuzzified, output.

3. DST output shows a good compatibility with well-test results in Burgan reservoir of "S" field, especially in wells 1, 3, 4 and 5. Therefore, well-test results of "S" field admit preciseness of DST output.

In "C" field, variation of testing results is not shown by DST output due to low variability of oil-rate values [0,650] 𝑏𝑏𝑙 𝑜𝑖𝑙 𝑑𝑎𝑦 .𝑚 in this field, comparing to "S" field, which is [0,4500]. As discussed before, cut-off method does not provide any relative pattern of productivity through pay zones.

4. DST method is free of cut-off determination, which is the heart of cut-off-based method of pay zone study. Cut-off determination is a challenging procedure, usually done by analogizing with nearby known oil fields or due to economic considerations. Therefore, proposed method gives a pure technical output, related to geological facts. From another viewpoint, generalizing DST method to adjacent fields or wells does not need another cut-off determination; so, easier generalization ability. In another word, no training is needed for the proposed method. Also there is a less need to optimization.

Conclusion

A novel procedure of net pay determination, based on Dempster-Shafer Theory (DST), is developed in this work. DST is a powerful decision-making tool in uncertain situations, like here in reservoir condition. Due to results, it is shown that DST-based method of net pay detection is capable in identifying productive zones through oil wells, which the output is compatible with output of conventional method and well-test data. The proposed methodology is admitted by cut-off method of pay determination in two reservoir types: carbonate (Mishrif) and clastic sandy reservoir (Burgan). In addition, by the means of the developed procedure, not only distinguished pay zones are detected fuzzily, which is very similar to real reservoir condition; but also the outputs are more compatible with well-testing results, comparing to conventional cut-off methodology. In addition, there is no need to train the proposed method. Even not an essential optimization is needed.

  information and interpretations, subsequently a big lack of knowledge in new areas. Even the same formation might show varying behavior in different locations or conditions! Therefore, in fact, "Net Pay Detection" is a general concept, containing several sub-problems: "Net Pay Detection: in Carbonates/ in Clastic reservoirs/ by using NMR log/ by Benefiting of Well Test Results/ … ".[START_REF] Worthington | The Role of Cut-offs in Integrated Reservoir Studies[END_REF] has listed 31 works on net pay determination, introducing studied reservoir rocks and contained fluid briefly. Set of reservoir rocks consists of ten types: sandstone, shaly sand, laminated shaly sands, laminated sandstone, silty sandstone, laminated muds and sands, tight sandstone, carbonate, limestone and dolomite. Set of fluid contents contains six types: oil, oil and gas, gas, gas condensate, heavy oil and extra-heavy oil. Surely, number of all possible conditions exceeds multiplication of ten by six, i.e. 60 situations! This number shows a big diversity in the problem of "net pay".

  Applying cut-off values on the well by rule of product Selecting one of water saturation cut-offs, which is more reasonable, as final cut-off of water saturation Determining separation points on histograms as potential cut-offs of water saturation Distinguishing different modes/distributions of water saturation by histograms Plotting histogram of water saturation in each well Determining borders of reservoir facies as cut-off values of porosity and kaolinite Distinguishing reservoir facies based on Kaol-Por plot Plotting Kaolinite-Porosity for each well

Figure 2

 2 Figure 2 and 3. Figure 2 shows net pay determination results in carbonate Mishrif reservoir of field "C";

Figure 2

 2 Figure 2 Net pay determination in carbonate Mishrif reservoir in seven wells of field "C". The vertical axe is depth

Figure 3

 3 Figure 3 Net pay determination in Sandy Burgan reservoir in five wells of field "S". The vertical axe is depth versus

Table 1

 1 List of articles, and introduced methodologies, datasets, instruments for net pay determination.

	Part	Article(s)	Methodology/ Datasets/ Instrument	Descriptions
	I:
	Part

Why to Solve a Single Problem by Multiple Solutions?

	Due to the literature, corresponding author have already proposed other solutions to the
	problem of "Net Pay Determination" based on different theories: Bayesian, fuzzy and diffusivity
	equation. There is a true critic: "Is it necessary to propose new methodology to solve a problem, which
	has already been solved by previous solutions?" In another word, wouldn't be redundant, approaching a
	single problem by multiple solution?
	Worthington, 2010)	Cut-off/ Shale volume, porosity and water saturation	Review paper on net pay
	(Mahbaz et al., 2011)	Optimization of cut-offs based on statistics, viscosity	Case study of Burgan reservoir in offshore
	(Masoudi et al., 2011)	Diffusivity equation/ Viscosity, compressibility, porosity and permeability from cores, conventional logs	Case study of Sarvak Formation, a carbonate reservoir
	(Masoudi et al., 2012a)	Bayesian classifier/ Shale volume, porosity, water saturation and well test	Case study of Sarvak Formation, a carbonate reservoir
			Case study of Sarvak
	(Masoudi et al., 2012c)	Fuzzy fusion/ LLD, LLD:LLS and well test	Formation, a carbonate
			reservoir
	(Masoudi et al., 2012b)	Comparing methods of cut-off, diffusivity, Bayesian and Fuzzy	Case study of Sarvak Formation, a carbonate reservoir
	Theory (DST), is developed in order to evaluate productive zones by belief function.

The majority of above mentioned articles are mainly based on cut-off method to determine net pays. There are only some novel methodologies in part III, which the bases are on data fusion approaches: Bayesian and fuzzy theories. The biggest pitfall of well-known cut-off procedure is that this method categorizes the whole gross into pay and non-pay zones that do not provide a realistic view of non-crisp nature of drilled well. Due to success of data fusion-based methods in assessing productive zones fuzzily in carbonates; in this study, another fusion-based approach, based on Dempster-Shafer

  Mishrif formation is overlaid on Khatiyah formation, and is distinguishable from overlaying Formation of Laffan by famous and distinct unconformity of Upper Cretaceous. The available dataset contains seven oil wells, drilled on "C" field, which is a relatively small oil field in the Persian Gulf; close to the border of Iran and Qatar. Zone of Oil-Water Contact (OWC) has been detected clearly in all the seven wells within Mishrif reservoir.

An important oil reservoir in the Persian Gulf and the Middle East, called Mishrif formation, equivalent to Upper Sarvak Formation with the age of Cenomanian-Early Turonian was available for this work. Based on sedimentologic studies, it is understood that Mishrif is deposited in a regressive cycle: starts with basinal (outer self) environment, followed by slope environment, continues with shoal or reefal environment and ends with lagoonal environment

[START_REF] Al-Dabbas | Depositional environments and porosity distribution in regressive limestone reservoirs of the mishrif formation, Southern Iraq[END_REF]

.In the studied oil field (name is not enclosed due to confidentiality, instead we call it Field "C" due to studying a carbonate reservoir through)

Acknowledgment

The authors would like to acknowledge Iranian Offshore Oil Company for permitting in publishing scientific achievements of the investigation. Also, thank to respected editor, Dr. Tahar Aïfa, and anonymous reviewers in spending time to review, and sharing their constructive ideas.