first break volume 22, November 2004
technical article
Quantifying uncertainty for mapping fracture intensity: an improved workflow
Patrick Wong and Sean Boerner*
Abstract Fractured reservoir characterization requires a good understanding of the spatial distribution of fracture intensity. In practice, the relationships between various geological, structural and seismic attributes and fracture intensity are complex and can be highly nonlinear. The complexity increases exponentially when there are many potential fracture-related attributes available, such as attributes derived from azimuthal AVO. This paper presents an improved integrated workflow for ranking the attributes, understanding data limitations, mapping fracture intensity with neural networks, and quantifying uncertainty. The study uses a set of soft attributes and cumulative gas production data, a proxy for fracture intensity, at the Pinedale Anticline at Wyoming to demonstrate the workflow and to show its usefulness in practice.
Introduction Modelling the spatial distribution of natural fractures in hydrocarbon bearing reservoirs can be a highly nonlinear problem. In such reservoirs, fluid-flow behaviour is governed by the geometries of the fractures (e.g. orientations, lengths and apertures), the intensity of the fractures for every fracture set and their connectivity. For flow simulation, it is also necessary to understand dynamic relationships of the fractures with the rock matrix (e.g. single/dual porosity, single/dual permeability) and apply a proper reservoir simulator for predicting the flow performance.
Fractured reservoir studies require the derivation of fracture intensity logs estimated from either static (e.g. image logs) or dynamic data (e.g. production data). As in the modelling workflow for conventional reservoirs, constraining the property models with soft data (seismic attributes) improves the predictions and lowers their uncertainty at regions with poor well control. The resulting fracture intensity map or volume can be used to constrain the spatial distribution of stochastic fracture planes, typically via the use of ‘discrete fracture network’ (DFN) modeling, and derive various flow properties (e.g. fracture permeability tensors, matrix-fracture interaction parameters) for each cell as inputs into a reservoir simulator (Ouenes and Hartley, 2000).
In this paper, we will focus on mapping the fracture intensity and quantifying the uncertainty of the predictions.
The improved workflow allows statistical ranking and selection of soft attributes and provides information on the inherent predictive power of the data set. We will use a standard back-propagation neural network to generate multiple stochastic maps of fracture intensity and then quantify the uncertainty in relation to the limitations of the data. We will demonstrate the workflow with an application to the Pinedale Anticline, Wyoming, USA.
Ranking attributes The use of regression-based mapping algorithms requires a training set with known input-output data pairs. For fracture intensity mapping, the output is the known fracture intensity data derived from the wells. The inputs are the soft attributes to be used to map the fracture intensity away from the wells. The types of attributes available vary from field to field. Some typical examples include azimuthal AVO attributes (e.g. crack density, Lynn et al. 1996), stress-related attributes (e.g. slopes and curvatures) and lithologies.
Ranking the relevance of the attributes is an essential step that allows us to select only the significant attributes for building the regression model. The reduction of input dimensions offers many advantages to numerical procedures, such as the increase of the number of training data per unknown in the model, the avoidance of over-constraining the model, and the reduction of CPU time.
In Wong and Boerner (2003b) we introduced a robust ranking algorithm known as ‘optimized piecewise rank correlation,’ or ‘OPRAC’ in short. This algorithm uses a rank correlation coefficient between each of the attributes and the output as a basis for ranking the significance of the attributes. It converts the data values to ranks, constructs optimal linear segments through the scatter-plot, and then calculates the rank correlation coefficient for each segment. All the coefficients are then aggregated and a weighted coefficient (ranging from 0 to 1) is obtained. Unlike the standard correlation coefficient approach that could give a very small coefficient for even highly correlated data (e.g. sine curve, symmetrical data), the proposed algorithm is able to capture any non-linear trends in the data. It also identifies dependent or cross-correlated attributes that add no additional information content to the regression model. This algorithm will be used in the paper.
Tel: 1-832-351-8372, Fax: 1-832-351-8776, Email: Patrick_Wong@veritasdgc.com
© 2004 EAGE
41
technical article
first break volume 22, November 2004
Figure 1 Indication of extrapolation of an input attribute (modified after Gauthier et al. 2000).
Attribute extrapolation While regression-based methods are powerful for integrating multi-attribute data, extrapolation occurs frequently in reservoir applications and cannot be ignored. In essence, there are two types of extrapolation. The first type is ‘spatial extrapolation’ that occurs when we need to make predictions away from the regions with good well control. This issue has been addressed extensively through the use of kriging variance that can be used to quantify the uncertainty of the predictions based on their separation distances from the hard data locations. The second type of extrapolation is ‘attribute extrapolation’ (Wong and Boerner, 2003a). This occurs when the full statistics of the input attributes are not present in the available training set to be used in regression.
The study presented by Gauthier et al. (2000) was one of the few published works that revealed the problem of attribute extrapolation. The authors examined the extent of extrapolation based on the coverage of each of the input attributes in the training set and the unsampled regions. Figure 1 shows a simple illustration of the limitation of an input attribute in the training set. As shown, the span of the training data does not cover the span of the attribute in the unsampled regions, and therefore extrapolation occurs for any regression-based methods. The authors also displayed a map of binary numbers indicating the areas of extrapolation if at least one out of the four main attributes experienced extrapolation. The resulting map showed that two-third of the entire field area experienced extrapolation. Note particularly that the indicator map is independent of the regression method chosen. Such a map serves as an excellent independent reality check and helps practitioners place confidence in the final results.
Coefficient of extrapolation In Wong and Boerner (2003a) we proposed the use of a ‘coefficient of extrapolation’ to fine-tune the binary extrapolation indicators as presented in Gauthier et al. (2000). We first examined the individual attributes and calculated the relative distance of departure from the extreme value (i.e. minimum or maximum) in the training data: the further away from the extreme value; the larger the coefficient. We then used the OPRAC coefficients to weigh the relative distances of all the attributes and an average value was calculated for each cell. This final value is the ‘coefficient of extrapolation’ ranged from 0 to 1. For the data within the coverage of the training data, the coefficient was set to zero (i.e. no extrapolation). Extrapolation is no longer a binary indicator, but it is now a matter of degree. Since the coefficient is available at each cell, the results can be plotted as a map or volume.
One deficiency of this method is that attribute extrapolation may occur in multidimensional space even though the individual attribute has no extrapolation problem. Figure 2 shows a cross-plot of two input attributes and two data points. Each attribute individually has no extrapolation for the two points, but extrapolation is clearly revealed when the data are posted on the two-dimensional input space. Therefore it is necessary to quantify the degree of extrapolation in the multidimensional input space.
In this paper we improved the calculation of the coefficient by searching for the closest training data point for each unsampled data point. The minimum distance is used to calculate the coefficient of extrapolation and was then normalized to 0 to 1. The procedure solves the problem presented in
Figure 2 Extrapolation in two input dimensions. The two circles represent two data points that do not experience extrapolation in the individual dimension, but they are outside the coverage of the training data in two dimensions.
42
© 2004 EAGE
first break volume 22, November 2004
technical article
Figure 2. If the unsampled data point were a training data point, then the calculated coefficient of extrapolation would be 0.0. If the coefficient were 1.0, the unsampled data point would be the furthest from all the available training data.
Mapping with neural networks Neural networks have become important in the geophysics community (Wong et al., 2002; Aminzadeh and de Groot, 2004). For fractured reservoir modelling, neural networks have been used with great success in various field studies (e.g. Zellou and Ouenes, 2001; Boerner et al. 2003). Their abilities to learn, adapt and generalize from data offer a number of practical advantages, including multi-attribute integration and uncertainty quantification with multiple stochastic realizations. For simplicity, this paper used a conventional backpropagation neural network.
In our previous works, the correlation coefficient between the actual fracture intensity values and the predicted fracture intensity values from the neural network was used to search for the favorable solutions. Our recent experience with a neural network application to a fractured reservoir study has prompted us to look for an alternative workflow. We found that if the available fracture intensity data was high skewed, it became difficult to find a satisfactory solution that conformed to our geological understanding of the field, even though the model converged in a stable fashion. The problem was that a few extreme values could skew the results and give a high correlation coefficient (Hirsche et al. 1997), producing a solution that would appear to be favourable, but would actually be poor if these extreme values were removed from the training set. These extreme values actually influenced the resulting solution much more than they should.
To circumvent the problem we improved our workflow by transforming the fracture intensity data into ranks. This essentially transformed the skewed distribution into a more uniform distribution, so all the training data points now have approximately an equal influence on the resulting correlation coefficient. The neural network modelling was performed in the rank space. When the model converged, the estimated ranks were back-transformed into the actual fracture intensity scale. With this rank transformation and back-transformation, we were able to obtain much more geologically sensible answers and therefore it is included in our improved workflow.
have been sparse until recently. The reservoir section is over 5000 ft (>1500 m) thick and composed of interbedded sands and shales deposited by a fluvial/alluvial system (Law and Johnson 1989). Matrix permeabilities are very low and in the order of microdarcies. Therefore, it is assumed that intersecting a connected fracture system is a requirement for profitable wells.
In this work we picked a top horizon with 19 wells within the area of interest and averaged the seismic attributes for the reservoir section. The four-month cumulative gas production data was used as a proxy for fracture intensity, because it is considered to be a good indicator of the expected ultimate recovery (Boerner et al., 2003). A total of 16 geological, structural (slopes and curvatures) and seismic attributes (attributes from AVAZ) were available as potential inputs into the neural network for mapping the fracture intensity.
Results and discussions
1) Ranking attributes The Pinedale data was first mapped on a 2D grid with 179 cells by 160 cells (28640 cells in total). The grid spacing was 200 ft (61 m) by 200 ft. The total areal coverage was over 41 square miles (106 km2). We used the OPRAC algorithm to rank the 16 attributes. The coefficients were ranged from 0.32 to 0.74. Since we had a small dataset, we selected only top four attributes for the modelling. Table 1 shows the relevance of these attributes with the OPRAC coefficients. The first, second, and third ranked attributes were the AVAZ attributes, and they were respectively the crack density, the coherency of the zero-incident shear wave (S-wave) seismic, and the azimuth of anisotropic gradient. The third ranked attribute was the top structure depth. Figure 3 shows the maps of the respective attributes and the well locations.
Table 1 The top four attributes with high correlation with fracture intensity.
Rank
1 2 3 4
Attribute
Crack density S-wave coherency Top structure depth Azimuth of anisotropic gradient
OPRAC Coefficient 0.74 0.58 0.50 0.44
Pinedale Anticline The Pinedale Anticline is located in Sublette County, southwestern Wyoming. A 3D seismic survey was acquired over the anticline and processed using Amplitude Versus Angle and Azimuth (AVAZ) analysis that estimates the local intensity and orientation of fracturing (Gray et al., 2003). These seismic attributes offer the advantage of dense lateral coverage and can be used to correlate with the well data. The area has been drilled for many decades, although profitable wells
2) Attribute extrapolation For the 2D map we compared the four attributes at the 28640 cells with the 19 training data and computed the coefficient of extrapolation map as shown in Figure 4. The red regions represent regions with significant extrapolation, while the blue regions have low extrapolation. Note that the coefficients were zero at the well locations. This figure indicated the inherent predictive power of the data
© 2004 EAGE
43
technical article
first break volume 22, November 2004
Figure 3 The top four attributes for fracture intensity prediction and the well locations for the study area of the Pinedale Anticline.
at hand. With this particular dataset, there are no regression-based methods that could reduce the extent of extrapolation. The only remedy is to modify the training set by: 1) using other attributes (probably less significant but with a wider coverage); and 2) adding synthetic or pseudo data to the training set.
Figure 5 shows a cross-plot of the first two top attributes (normalized to standard scores, or Z-scores), crack density and S-wave coherency, for both training and unsampled data. As displayed, the training data covers only a small part
of the overall data coverage. We also show a cross-plot of the resulting coefficient of extrapolation against the S-wave coherency in Figure 6. The S-wave coherency with data greater than 2 experienced significant extrapolation.
3) Mapping fracture intensity We used the four top attributes as inputs to the neural network models. All the 19 fracture intensity values from the wells were transformed into ranks. Ten models were built with different initialization settings and resulted in 10 differ-
44
© 2004 EAGE
first break volume 22, November 2004
technical article
Figure 4 A map of coefficient of extrapolation.
Figure 6 A cross-plot of coefficient of extrapolation and Swave coherency.
Figure 5 A cross-plot of S-wave coherency and crack densi-
ty for training and unsampled data.
Figure 7 Average fracture intensity map.
ent stochastic realizations. The convergence was fast and stable. The results were then back-transformed into fracture intensity. We calculated the corresponding average and standard deviation (STD) maps from the 10 realizations. The results are shown in Figures 7 and 8 respectively. Note that the regions with high fracture intensity followed the same patterns of the crack density map as shown in Figure 3a, because it has the highest OPRAC coefficient. In cases where the natural fractures are more or less aligned in the same direction, the crack density from AVAZ can be used as an excellent proxy for estimated fracture density.
4) Uncertainty quantification The low STD regions, shown as blue in Figure 8, represent little variations between the realizations. These regions are often treated as low risk areas. What we must be cautious about is that the low STD areas represent the areas in which the regression algorithm itself is confident in the predictions, not the degree of extrapolation from the training data that was used to generate the predictions. Therefore it is important to conduct a reality check with the computed uncertainty as shown in the STD map.
One major use of the coefficient of extrapolation map is to crosscheck the uncertainty estimates from stochastic sim-
© 2004 EAGE
45
technical article
first break volume 22, November 2004
Figure 8 Fracture intensity standard deviation map.
Figure 10 Fracture intensity standard deviation map with coefficient of extrapolation greater than 0.11.
polation we used a coefficient of extrapolation of 0.11 (i.e. the highest 20%) as a cutoff value and displayed the STD map with coefficients greater than 0.11 in Figure 10. As shown, for these areas with significant extrapolation, the resulting high STD may be deemed as reasonable predictions while we may remain sceptical about the predictions at low STD. Therefore the coefficient of extrapolation map may be used to cross-check and validate prediction variability from the stochastic simulations and improve our ability to quantify uncertainty. These maps may be used to locate potential pseudo wells in order to obtain data for improving the information content and thus bolstering the statistics of the training set. They may further improve infill-drilling decisions and indicate potential sweet-spots with higher confidence.
Figure 9 A cross-plot of fracture intensity standard deviation with coefficient of extrapolation.
ulation. Figure 9 shows a cross-plot of fracture intensity standard deviation with the coefficient of extrapolation. This plot reveals that some small STD values (e.g. <20000) were in fact derived from data with very low predictive power (i.e. high coefficient of extrapolation). Therefore there is significant extrapolation in the data, but neural networks consistently give similar predictions.
In order to show the regions of high coefficient of extra-
Conclusions In this paper, we introduced an improved workflow for mapping fracture intensity from well data. The workflow started with use of an improved algorithm, ‘optimized piecewise rank correlation’ (OPRAC), for ranking the relevance of the soft attributes and then mapping fracture intensity with neural networks in the rank space. ‘Coefficient of extrapolation’ was proposed to understand the inherent predictive power of the data at hand and to crosscheck and validate the variability of the predictions from neural networks. The workflow was demonstrated with the data from the Pinedale Anticline in Wyoming. The study shows the usefulness of the workflow and demonstrates how it can be applied in practice.
Acknowledgements The authors would like to thank the anonymous referees for
46
© 2004 EAGE
first break volume 22, November 2004
technical article
reviewing this work and Veritas DGC for permission to publish this work.
References Aminzadeh, F. and de Groot, P. [2004] Soft computing for qualitative and quantitative seismic object and reservoir property prediction. Part 1: Neural network applications. First Break 22, March, 49-54. Boerner, S., Gray, D., Zellou, A., Todorovic-Marinic, D. and Schnerk, G. [2003] Employing neural networks to integrate seismic and other data for the prediction of fractures intensity. 2003 SPE Annual Technical Conference and Exhibition, Denver. SPE 84453. Gauthier, B.D.M., Zellou, A.M., Toublanc, A., Garcia, M. and Daniel, J.-M. [2000] Integrated fracture reservoir characterization: A case study in a North Africa field. SPE European Petroleum Conference, Paris, SPE 65118. Gray, F.D., Todorovic-Marinic, D. and Lahr, M. [2003] Seismic fracture analysis on the Pinedale Anticline: Implications for improving drilling success. EAGE 65th Conference and Exhibition, Stavanger, paper C-20. Hirsche, K., Boerner, S., Kalkomey, C., and Gastaldi, C. [1998] Avoiding pitfalls in geostatistical reservoir characterization: A survival guide. The Leading Edge, April, 493-504. Law, B.E. and Johnson, R.E. [1989] Structural and strati-
graphic framework of the Pinedale Anticline, Wyoming, and the multiwell experiment site, Colorado. In: Geology of Tight Gas Reservoirs in the Pinedale Anticline area, Wyoming, and the Multiwell Experiment Site, Colorado, (ed. B.E. Law and C.W. Spencer), B1-B11, U.S. Geological Survey Bulletin 1886. Lynn, H.B., Simon, K.M. and Bates, C.R. [1996] Correlation between P-wave AVOA and S-wave traveltime anisotropy in a naturally fractured gas reservoir. The Leading Edge 15, 931-935. Ouenes, A. and Hartley, L.J. [2000] Integrated fracture reservoir modeling using both discrete and continuum approaches. SPE Annual Technical Conference and Exhibition, Dallas, TX, SPE 62939. Wong, P., Nikravesh, M. and Aminzadeh, F. [2002] Soft Computing for Reservoir Characterization and Modeling, Physica-Verlag, Heidelberg, 586 pp. Wong P. and Boerner S. [2003a] Fracture intensity modeling using soft computing approaches. SEG Summer Research Workshop on “Quantifying Uncertainty in Reservoir Properties Prediction,” Galveston, TX, abstract only. Wong, P. and Boerner, S. [2003b] Ranking geological drivers: A comparison study. Computers & Geosciences 30, 91-100. Zellou, A. and Ouenes, A. [2001] Integrated fractured reservoir characterization using neural networks and fuzzy logic: Three case studies. Journal of Petroleum Geology 24, 459-476.
EAGE Awards 2004
EAGE
MORE DETAILS?
WWW.EAGE.ORG
Every year EAGE presents several prestigious awards to members of the Association who have made a highly significant contribution to one or more disciplines, to authors of the best papers published in EAGE journals and also to presenters of the best oral and poster presentation at the Annual EAGE Conference.
• Do you know some outstanding contributors to geoscience? • Have you read some exceptionally high standard articles this year?
Then make a nomination for one of the EAGE Awards!
Please send your nominations before the end of November 2004 to eage@eage.org
© 2004 EAGE
47