Chapter 7 Mapping geological properties and uncertainties

Chapter 7 Mapping geological properties and uncertainties

CHAPTER 7 Mapping Geological Properties and Uncertainties C O M P U T E R CONTOURING Some of the geological properties that we use to define prospec...

2MB Sizes 0 Downloads 88 Views

CHAPTER

7 Mapping Geological Properties and Uncertainties

C O M P U T E R CONTOURING Some of the geological properties that we use to define prospects are continuous surfaces, such as the contact between two stratigraphic units (expressed as either the top or bottom of an interval) and the contacts between oil and water and gas and oil in a reservoir. We may know the elevations of formation tops in wells that penetrate the formations, or we may have estimates of their depths at locations along a seismic survey. However, we obviously do not know the true values of such surfaces at a prospect site because this would mean that the prospect has already been drilled! We must, in some manner, calculate estimates of the elevations of the surfaces at the prospect locality. More generally, we would like to estimate values of the surfaces of interest at all locations within the region, because we could then represent the surfaces as contour maps. The construction of contour maps was one of the first major tasks to be delegated to computers by exploration geologists; primitive contouring programs appeared in the early 1960's (c/. IBM, 1965). Initially, programmers were happy if their code could generate an acceptable-appearing contour map in a tolerable length of time. By present standards, the mainframe computers available in the 1960's and 70's were enormously expensive, had minuscule memories, and were agonizingly slow. Only the larger oil companies, universities, and government research laboratories could aff^ord to experiment with computer contouring. Because of the limited capabilities

Computing Risk for Oil Prospects — Chapter 7 of the available machines, great emphasis was placed on computational efficiency and the use of programming "tricks" that would allow larger data sets to be processed than would otherwise be possible. The major oil companies realized that by automating one of the geologist's most time-consuming tasks, computer contouring could greatly speed the work of their exploration groups. In a competitive business like oil exploration, being first is often critical. In addition, these same companies had already invested heavily in computers for geophysical processing, record-keeping, and business purposes, and it seemed natural that their explorationists should make use of this investment as well. The geologists, however, often took a dim view of computer contouring and considerable effort was devoted to convincing the skeptics that computer-generated maps were acceptable (Dahlberg, 1975; Walters, 1969). Commercial software houses developed contouring programs to appeal to the large oil companies with their enormous investments in computing equipment. Like the computers themselves, the programs were large, complex, and extremely expensive. Because of their complexity, the programs were not operated by geologists, but by computer specialists. Since the programs were so expensive to purchase and costly to maintain and run, it was essential that the visual quality of their output reflect the magnitude of the investment. Consequently, great emphasis was placed on graphical perfection and cartographic embellishments. Also, because the programs were used in conjunction with large corporate and commercial well data bases, they included sophisticated data management capabilities. The one critical aspect of contouring that was overlooked during this period of software development was an assessment of the reliability of the finished computer-drawn contour maps. Maps produced by computer were judged on the basis of their visual appearance, which meant how closely their features resembled those on manually contoured maps. About the only quantitative measure of acceptability was the degree to which a map "honored the data points." Arguments that this was not a valid criterion for judging a mapping procedure whose primary function was estimation and not reproduction were roundly ignored (Davis, 1976). There are extensive discussions about how computer contouring programs operate, and exhaustive comparisons of the characteristics of different algorithms. We will not pursue these details, but instead refer the reader to discussions in Davis (1986), Hamilton and Jones (1992), Jones, Hamilton, and Johnson (1986), Robinson (1982), and especially Watson (1992), who has a complete bibliography on contouring algorithms. Here we will briefly outline the underlying principles of computer contouring and place them in the broader scope of geostatistics, which provides a means of assessing the reliability of contour maps in probabilistic terms. 138

Mapping Properties and Uncertainties In risk assessment, it is essential that we regard contouring as more than a quick and convenient way of creating a pretty picture of, say, the subsurface structural form of a prospect. We must be able to use contouring as a forecasting tool, producing an estimate of the likely shape of a structure that is only dimly perceived from scattered observations. Furthermore, we must be able to assess the uncertainty in our estimate, and to define the structure's likely upper and lower limits. Only when we do this can we include the potential range of the prospect in our appraisal.

HOW CONTOUR MAPS ARE MADE All computer contouring algorithms make several assumptions about the surface being mapped. The surface is presumed to be single-valued at each point or geographic location, to be continuous everywhere within the limits of the map, and to be autocorrelated over a distance that is greater than the typical distance between the available data points. ("Autocorrelation" is a statistical concept borrowed from time series analysis; it indicates the degree of similarity between a signal and itself after it has been shifted.) If these assumptions are valid, the known values of the surface at the control points can be combined to produce estimates of the surface between the control points. A geological surface such as the contact between two stratigraphic intervals obviously is single-valued and continuous unless it has been thrust faulted or recumbently folded. Other geological properties of interest, such as porosity, are not so obviously single-valued since more than one measurement of porosity may be available for a reservoir interval in an individual well. However, we can think of such multiple observations as samples that can be represented by a statistical value such as their mean. We may assume that the mean porosity of an interval is a single-valued property that is continuous, and hence mappable. Although geophysicists routinely utilize the concept of autocorrelation in seismic processing, it is less familiar to geologists. However, the concept is essential for an understanding of spatial phenomena. "Spatial autocorrelation" refers to the degree to which values at specific locations are related to values at other locations which are separated by a constant distance in a fixed direction. It provides a formalization of the commonsense concept that points on a surface are very much like nearby points, and less like more distant points.

Conventional Contouring Programs Computer contouring programs perform a number of distinct operations. First, the program must sort through the data and organize them so data 139

Computing Risk for Oil Prospects — Chapter 7 points can be selected quickly according to their location. Next, the program establishes an imaginary regular gridwork of locations across the map area where the estimates of the surface will to be made in succeeding operations. Ordinarily, the user must specify the dimensions of this grid, either as the number of rows and columns or the distances between the rows and columns. Other information may be necessary if the map is to be either larger or smaller than the extremes of the locations of the data points. Once this gridwork is defined, the program is ready to estimate the values of the surface at the nodes of the grid. To estimate an individual grid node, the nearest control points—usually well locations or seismic shot points—are combined in the form of a weighted average. At this time, the user may be faced with a number of decisions. How is "nearest" to be defined? How many "nearby" points should be used in each estimate? And what weights should be applied to each point? Diff^erent programs incorporate different choices in their coding, and the vendors of some commercial software packages tout their particular choices as superior. Other programs leave these decisions up to the user, and provide a variety of diff'erent alternatives. The fact is, such choices are arbitrary; the combination of data selection criterion, number of points, and weighting function that proves eff^ective with one data set may not be best for a diff^erent variable or a diff'erent set of data points. The algorithms are sensitive to the characteristics of the property being mapped and the geometric arrangement of the observations. Experience has shown that some algorithms tend to work well with the type of data that commonly are mapped in petroleum exploration, such as the structural configuration of a subsurface horizon or the thickness of an interval. Typically, the program will select 8 or 16 control points around a grid node that is to be estimated. The area around the grid node is divided into eight wedge-shaped sectors and the nearest one or two points are found in each sector. This ensures some degree of radial control of the estimate even if the data points are closely spaced along widely separated lines (as is typical of seismic and airborne geophysical measurements). Finally, the selected points are combined as an average in which each observation is weighted by the inverse of the square of the distance from the observation to the grid node. The weights are scaled so the most distant point has a weight near zero and the sum of the weights is equal to 1.0. There are many other factors that must be considered in practice. A maximum search radius and a minimum acceptable number of control points for each estimate must be specified, or otherwise the algorithm may incorporate far distant control points in areas where there are few observations, or fail to estimate grid nodes near the edges of the map. The literature on contour mapping (which is extensive but widely scattered) is 140

Mapping Properties and Uncertainties Table 7.1. Summary statistics for dry holes and producing wells in the Magyarstan training area. Data are contained in the file TRAINWEL.DAT on the diskette. Elevation Thickness Shale Bedding ratio index (m) (m) Producing wells (n = 18) Maximum -1404.0 90% -1409.4 -1431.2 75% 50% -1554.5 25% -1676.0 10% -1682.7 -1689.0 Minimum Mean -1550.9 Standard deviation 119.2

46.73 46.52 44.77 38.97 36.57 35.28 35.09 40.73 4.29

0.408 0.390 0.379 0.371 0.342 0.315 0.293 0.361 0.029

0.270 0.260 0.250 0.243 0.235 0.225 0.215 0.243 0.013

Dry holes (n == 65) Maximum 90% 75% 50% 25% 10% Minimum Mean Standard deviation

46.77 42.52 40.69 37.72 34.28 32.17 27.61 37.26 4.19

0.650 0.567 0.514 0.476 0.445 0.407 0.392 0.483 0.056

0.338 0.319 0.295 0.275 0.243 0.231 0.196 0.274 0.033

-1190.0 -1355.4 -1468.5 -1557.0 -1618.5 -1659.0 -1705.0 -1529.6 114.4

filled with discussions of the effects that these and other design decisions have on the resulting maps. Table 7.1 gives statistical summaries of measurements on four geological properties of the XVa Limestone as measured in 83 wells drilled in a small area of Magyarstan on the Zhardzhou Shelf that we refer to in this book as the "training" area. The data are included in the diskette file TRAINWEL.DAT. The geological setting and significance of these particular variables are described in Chapter 1, and Bayesian estimates of conditional probabilities based on the magnitudes of some of these variables are given in Chapter 5. Here, we will use these properties to illustrate computer contouring and subsurface analysis performed using the RISKMAP software, which is included on the accompanying diskettes. 141

Computing Risk for Oil Prospects — Chapter 7 The area to be mapped measures 35 x 35 km. Each geological property will be represented by a grid of values that contains 36 rows and 36 columns, so the contouring program must estimate values of a surface at grid nodes spaced at 1-km intervals. The origin of the grids is in the southwest corner and has coordinates of 1740 (east-west or Xi-direction) and 6455 (north-south or X2-direction). These are transformations of realworld UTM coordinates. We will use a simple contouring algorithm that estimates each grid node from the eight nearest wells, weighting the value at each individual well by the inverse of the square of the distance between the grid node and the well. Figure 7.1 is a contour map of the subsurface structural configuration of the top of the XVa Limestone unit; the contours represent meters below sea level. Contour lines have been drawn at intervals of 25 m. In addition to local irregularities in form, there is a general tendency for the surface to dip to the southeast. Figure 7.2 is an isopach map showing the thickness of the XVa Limestone unit, measured in meters. Thicker intervals are interpreted as being carbonate buildups or reefs, particularly if the carbonate interval is massive (see Chapter 1). Figure 7.3 is a contour map of the shale ratio, based on the average response of the gamma-ray log. Lower values represent relatively pure limestone that is free of included shale. Figure 7.4 is a contour map of the bedding index, a measure in which low values represent massively bedded limestones and thick shales at one extreme and high values represent finely interbedded limestones and shales at the other extreme. The bedding index is derived from the standard deviation of the gamma-ray log, which is sensitive to the shale content of the rocks.

Trends and Residuals Sometimes it is helpful to transform a contour map in different ways, to clarify features or to separate the map into components in order to see relationships that might otherwise be hidden. A widely used procedure is trend surface analysis, in which a dipping plane or simple curved surface is used to approximate the observations. This "trend" is fitted by a least squares criterion, in which the coefficients of the surface are determined in a manner that minimizes the sum of the squared deviations of the actual observations from the trend surface. This is exactly the same procedure that is used to fit a line to data in a cross plot, as described in Chapter 4. The fitted trend is an approximation of the data, and is a smooth representation of the large-scale form of the surface. If the observations are structural elevations, the trend often is interpreted as a regional component of the structure. Davis (1986) provides an extensive discussion of trend surface analysis and its alternative interpretations. 142

Mapping Properties and Uncertainties

6490.

6485.

6480. K

6475.

6470. r

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.1. Structure contour map of top of XVa Limestone in the Magyarstan training area. Coordinates given in kilometers. Contours given in meters below sea level. Contour interval is 25 m. Solid dots are producing wells; other symbols indicate dry holes. Summary statistics for the 83 wells are listed in Table 7.1.

The trend can be subtracted from the actual surface, leaving residuals where the two do not coincide. These residuals also can be displayed in the form of contour maps, consisting of positive "highs" where the actual surface lies above the trend, and negative "lows" where it is below the trend. In effect, the residuals represent small-scale variations in the shape of the mapped property that are not contained in the large-scale trend. If we regard the trend as representing aregional component of structure, the residuals represent local features. 143

Computing Risk for Oil Prospects — Chapter 7 6490.

6485.

6480. h

6475.

6470.

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.2. Isopach map of the thickness of the XVa Limestone in the Magyarstan training area. Contour interval is 2 m.

A trend surface is calculated from the original observations using a statistical procedure that is a direct extension of the regression technique discussed in Chapter 4 and incorporated in RISKSTAT. The mapped variable, perhaps elevations of a subsurface horizon as measured in a number of wells, is the dependent Y variable and the geographic coordinates of the well locations are two independent X variables. A linear trend is defined by the equation Y = a -\- b^Xi -\-b2X2, where Y might be elevation, Xi the east-west coordinate, and X^ the north-south coordinate. The coefficient 61 gives the slope of the trend in the east-west direction, and 62 gives its slope in the north-south direction. The fitted surface has the form of a uniformly dipping plane. Coefficient a is the intercept, or value of the 144

Mapping Properties and Uncertainties

6490.

6485.

6480.

6475. r

6470. h

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.3. Contour map of shale ratio, calculated as the average gamma-ray response in the XVa Limestone in the Magyarstan training area. Contour interval is 0.02 units. Low values indicate intervals with little clay content, high values indicate intervals with relatively high clay content. trend at the origin, and reflects the average value of y . The coefficients are found by solving a set of simultaneous equations that are direct extensions of those described in the section on fitting lines in Chapter 4. If a dipping plane seems too simple a description of the regional trend, we can expand the trend surface equations to fit more complicated surfaces. This is done by adding new variables to the trend surface equation and solving for the additional coefficients by expanding the set of simultaneous equations. The new variables are powers and cross products of the geographic coordinates. For example, a second-degree trend surface has an equation: Y = a-[- hiXi + 62X2 + h^Xl -f- 64XI + 65X1X2. The additional 145

Computing Risk for Oil Prospects — Chapter 7 6490.

6485.

6480.

6475.

6470. h

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.4. Contour map of bedding index calculated as standard deviation of gamma-ray response in the XVa Limestone in the Magyarstan training area. Contour interval is 0.02 units. Low values indicate massive intervals, high values indicate thin-bedded intervals. squared terms allow the surface to bend or change slope in each of the two coordinate directions. So, a second-degree trend surface might have the form of a dome, basin, or saddle. More complicated forms are possible by using additional terms of higher powers, such as Xf, X^, or combinations of powers of the two coordinates, such as XJX2' Once the trend surface coefficients have been estimated from the data, the trend equation can be quickly evaluated for any set of geographic coordinates. Therefore, it is a simple matter to evaluate the equation at all the grid node locations where our contouring program has made an estimate of the surface. We will then have two grids of values, one containing estimates 146

Mapping Properties and Uncertainties 6490.

6485.

6480.

6475.

6470.

6465. r

6460.

6455.

1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.5. Second-degree trend surface fitted to elevation of top of the XVa Limestone in 83 wells in the Magyarstan training area. Contours given in meters below sea level. Contour interval is 25 m.

of the surface itself, and the other containing the trend. By subtracting the two grids, we produce a grid of residuals. All three of these grids can be displayed as contour maps. Figure 7.5 is a second-degree trend surface of the subsea elevation of the top of the XVa unit, the surface shown in Figure 7.1. The pronounced structural dip to the southeast is readily apparent. Figure 7.6 is a map of residuals produced by subtracting the grid of Figure 7.5 from that of Figure 7.1. The regional dip has been removed, leaving local areas that are relatively "high" with respect to the trend distinguished from areas that are relatively "low." 147

Computing Risk for Oil Prospects — Chapter 7 6490.

6485.

6480.

6475.

6470.

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.6. Residuals from second-degree trend surface (Fig. 7.5) fitted to elevation of top of the XVa Limestone in the Magyarstan training area. Residuals are found by subtracting trend surface grid from structural grid used to construct Figure 7.1. Contour interval is 10 m. Often a structural contour map, such as Figure 7.1, will show no traces of structural closure, although "noses" and other plunging features may be evident as deflections in the contour lines. If the regional trend is removed, these deflections may be transformed into closed features. These closed contours are easier to discern than are the deflections from which they are derived, and in some circumstances they even can be interpreted as originally closed features that have been tilted by regional movement. However, it's not necessary to attach a genetic significance to a trend surface and its residuals; it is sufficient to regard them as arbitrarily transformed versions of the original map that may be useful in the search for oil and gas. 148

Mapping Properties and Uncertainties Trend and residual maps can be prepared for any geological variable, not just elevations of subsurface stratigraphic units. Residuals from an isopach map show locations where the mapped unit is unusually thick or thin, and trend maps of porosity and other highly erratic properties may be easier to interpret than contour maps of the raw data. Just as there are many geologic properties that could be mapped, there are any number of trend surface and residual maps that could be made. The critical issue is whether such maps are useful in delineating prospective areas.

GEOSTATISTICS IN RISK ASSESSMENT Conventional contouring programs may produce acceptable maps of a geological surface, in the sense that the appearance of the maps meets our expectations. Unfortunately, there is no way to test the overall "goodness" of a conventionally produced contour map, nor is there any way to delineate areas on a map that may be more uncertain or less reliable than other areas. However, it seems reasonable that predictions about a subsurface unit should be better where the unit is penetrated by many wells than where it is relatively undrilled. Because conventional contouring algorithms are empirically based, there is no underlying theory that can be used to produce measures of reliability or uncertainty. Geostatistics is a branch of applied statistics that treats the variation of properties through space, in both two and three dimensions. As you might expect, the concept of autocorrelation plays a key role in this field. Geostatisticians have developed methods for estimating the spatial autocorrelation of a surface from scattered observations such as exploratory drill holes, and then using models of the autocorrelation to construct optimal estimates of the surface. These estimates can be displayed in the form of contour maps. Conventional contouring programs calculate estimates of a surface as distance-weighted averages of nearby points. The weights are assigned arbitrarily, usually as the inverse of the squared distances between the location where the estimate is being made and the control points being used. In contrast, geostatistical estimators are "custom-made" for the degree of autocorrelation in the actual surface, and the specific arrangement of points around every location being estimated. The geostatistical literature is vast and often couched in nomenclature and mathematics of fearsome complexity; Isaaks and Srivastava (1989) provide a modern, comparatively lucid introduction to the topic, while Cressie (1991) examines the field from the viewpoint of a statistician. A classical treatment with applications to mining is given by Journel and Huijbregts (1978). Deutsch and Journel (1992) provide software. The seminal works by Matheron (1962, 1965) are reserved for the mathematically masochistic. 149

Computing Risk for Oil Prospects — Chapter 7

T h e Semivariance The spatial autocorrelation of a surface is measured by a special statistic called the semivariance, which is simply the average squared difference between pairs of points that are separated by a constant distance. This statistic can be calculated for various distances and the results plotted as a graph of distance versus semivariance. Such a graph is called a semivariogram. There should be no difference between a measurement made on a surface at a point and a second measurement made at the same point (that is, between two measurements that are separated by a distance of zero), so the semivariance for a surface over a distance of zero should be zero. In other words, the semivariogram plot should go through the origin. Since values at locations that are separated by small distances should be similar, their average difference will be small and the semivariance of a surface over short distances should be low. If more distant locations on the surface are compared, there will be greater differences between values on the surface, and the semivariance over large distances will be greater. However, the differences between pairs of points will not continually increase with distance, but instead will become a more or less constant value. This upper limit will be numerically equal to the variance of all measurements made on the surface, without considering their spatial locations. A simple equation for semivariance is

^^^

2^

where xi is the value on a surface at some point i and xi^h is the value at another point located a distance h away. There are n such points, so there are 2n possible pairs of comparisons. The semivariance corresponding to a distance h between the pairs of points is indicated by 7/1. However, it is obvious that calculating numerical values for 7/^ requires that we have a set of points on the surface that are separated by a constant distance. If we have such a traverse of equally spaced points (such as a seismic hne or row of development wells in a large field), we can easily calculate the semivariance for multiples of the spacing between observations (Fig. 7.7). The multiples of spacing are referred to as lags, a terminology borrowed from time series analysis. Often, the X axis of the semivariogram is given in terms of lags, or multiples of the basic distance between pairs of points, rather than in actual distances. Figure 7.8 is a plot of average porosities measured in a series of production wells drilled into a shallow sandstone reservoir in eastern Kansas. The field is drilled with a 10-acre spacing, so the 20-well traverse extends for 2.5 mi and the unit distance between wells is about 660 ft. The first step 150

Mapping Properties and Uncertainties

h= 1

+ + + + + + + + + + + + + + + + + +

\J\JKJKJ\JKJ\J\J h=2

+ + + + + + + + + + + + + + + + + +

Figure 7.7. Straight traverse along a row of equally spaced wells showing pairwise comparisons between wells. For h = 1, every well is compared with its neighbors; for h = 2, every well is compared to every other well; and for h = 3, comparisons are between wells separated by two intervening wells. in calculating the semivariance is to find the squared differences between all possible pairs of points for successive lags. These differences are shown in Figure 7.9, plotted against lags. The actual distances can be found by multiplying the lag number by 660 ft, the basic distance between wells. (Although there appear to be fewer measurements at small lags, this is because the plotted points overlap.) Squared differences are shown for lags from 1 through 10, a distance of 1 mi. Next, averages are found for each of the lag distances. In Figure 7.10, the averages are shown by heavy lines and "box-and-whisker" plots show the distributions of the individual squared differences for each lag. It's apparent that both the averages and the scatter in the squared differences increase up to about lag 6 and then remain more or less constant. Figure 7.11 shows the completed semivariogram, without the intermediate values used in its calculation. This is referred to as an experimental semivariogram because it is based on a sample of observations rather than on theoretical considerations. The semivariogram will be used for calculation of weights in the geostatistical contouring algorithm, which goes by the name of kriging. (The name honors a prominent South African mining engineer, Danie Krige, who introduced statistical methods into mine evaluation.) In kriging, the semivariance may be required for any distance, not just those corresponding to the discrete lags of the experimental semivariogram, so we must create a continuous model of the semivariogram. This model semivariogram is an 151

Computing Risk for Oil Prospects — Chapter 7 25-

I

I

I

I

10 Observation

I

I

I

1

15

1

I

20

Figure 7.8. Profile showing clianges in average porosity along a traverse through 20 wells in an eastern Kansas oil field.

150-

§100-

i •8 (0 CO

-

J! 1

L '2

'3

' 4

' 5

' 6 Lag

' 7

' 8

' 9

• 10

Figure 7.9. Squared differences between pairs of wells separated by 1 to 10 lags along traverse shown in Figure 7.8. Lag interval is 660 ft. 152

Mopping Properties and Uncertainties IOU~ •

__ •

§1000)

1

1



J- 50-

-h





1

1

1

1

1

1

rn

• O"'

f

-1

1

m 1

n 11

11

r

r

JJJ

H

LfJ

ffl

m m m

' :2^~ 3 ' -4 ^ 5 ^ (5

iTi

r 11ij

1

il

bb 1

h"

.r

+ 1

MT—r

1

3

1 1

-

7 '



B '

9 ' 10 '

Lag Figure 7.10. Distributions of squared differences between pairs of wells shown by "box-and-whisker" plots. For each lag, "boxes" enclose central 75% of squared differences, and "whiskers" extend to cover 90% of squared differences. Line connects average squared differences, which are equivalent to the semivariance. Lag interval is 660 ft. idealized representation of the spatial continuity of the surface. Because the model has the form of an equation, it can be evaluated for any required distance. Geostatisticians have developed a library of model equations that are especially useful in representing the experimental semivariograms that are found commonly in petroleum exploration. Some of these are shown in Figure 7.12; all begin at the origin and rise smoothly to an upper limit called the sill. They reach this limit at a distance called the range. The values at points on the surface which are closer to each other than the range are related; if the points are separated by distances greater than the range they are statistically independent of one another. The experimental semivariogram may not be so well-behaved as the one shown in Figure 7.11, especially if the number of wells is hmited. The semivariance may change erratically at successive lags, making the task of fitting a model seem hopeless. In such circumstances, the best that can 153

Computing Risk for Oil Prospects — Chapter 7

Figure 7.11. Experimental semivariogram for average porosity along traverse through eastern Kansas oil field. Lag interval is 660 ft.

be done is to define a model that approximately forms a convex hull or envelope around most of the semivariances, and which reflects as closely as possible those nearest the origin. Semivariogram models cannot be fitted readily by least squares or similar techniques, because the objective is not to find a model that goes through the middle of the scatter of experimental values, but rather one that encloses most of the values. This will result in a conservative estimate of spatial continuity and conservative estimates of the amount of uncertainty about the surface. The simplest model for the semivariogram is a straight line, called a linear model (Fig. 7.12a). It begins at the origin and rises at a constant rate until it reaches the sill; from that point it is a constant. In equation form, 7/i = ah

forh
Ih = ^ 0

{ox h> a

where CTQ is the variance of all of the points measured on the surface. The model is fitted by trial and error, adjusting the slope, a, of a straight line 154

Mapping Properties and Uncertainties u

\ - < - '

O

o

o

^^^^ "^"^O

a D^

\ - < -

a

o

^--^-"^

\j

o

-'"O

0 O

a B^

c

b

V-

CO

> *E o

L..- ^2 .... "-a— —<^

%

CO

6"

o O

a D^

L...^....

c o

o

o

" ^ . y ^ ^

a »^

d

Distance, h Figure 7.12. Popular models for semivariance (theoretical semivariograms). Range is indicated by a, value of sill by G\. (a) Linear semivariogram with a constant slope up the sill, (b) Spherical semivariogram. (c) Exponential semivariogram. (d) Gaussian semivariogram. through the origin until it just encloses the majority of the semivariances that are near the origin. The most widely used model for the semivariogram is the spherical model (Fig. 7.12b) which is very straight near the origin but bends and merges smoothly with the horizontal line that represents the sill. The equation for the spherical model is 2/3/1

h^\

7/i = ^0 loxh> a To fit a spherical model, all that is necessary is to define a range, o, for the semivariogram, and then to calculate 7/^ for various values of h. As a first 155

Computing Risk for Oil Prospects — Chapter 7 approximation, the sill may be set equal to the variance of the observations. It may be necessary to experiment with different levels for the sill and values of a in order to find a semivariogram that is acceptable. Sometimes the experimental semivariogram rises faster at the origin than does the spherical model. In such instances, an exponential model (Fig. 7.12c) may provide a better fit. The equation for this model is 7;, =al{l-

e"^/^)

7/1 = o^l

for / i < a for h > a

Again, only the range, a, and the sill, aQ, are needed to fit this model. The range has a slightly different meaning in this model, because the exponential curve approaches the value of the sill asymptotically and so never quite becomes equal. To avoid this problem, the range is defined based on a nominal distance at which the semivariance "approaches closely" to the sill. A practical rule is to set the range as y/Sa, where a is the nominal distance at which the semivariance reaches 95% of the value of the sill. Experimentation with different sills and values of a may be needed to fit the model to the data (Olea, 1991). Geologic properties that are exceptionally smooth and continuous will produce an experimental semivariogram that rises slowly at low lags, then more rapidly with increasing distance. Such behavior can be represented by a Gaussian model (Fig. 7.12d), which is parabolic in form near the origin. The Gaussian model equation is 7/1

= al{l-

7/i =

(JQ

e-^'^'/a^)

for / i < a for h>

a

Like the exponential semivariogram, the Gaussian model approaches the sill asymptotically. The same procedure can be used to select an appropriate value for the range. Other more complicated models are found in the geostatistical literature, as well as models formed by combining simpler forms; these we leave to the geostatistical sophisticate. We will only mention that sometimes the experimental semivariogram does not seem to go through the origin, either because of fluctuations in the surface over distances shorter than the spacing between wells, or because repeated measurements at a location may result in different values (properties such as porosity may exhibit this behavior). This is referred to as a nugget effect, and simply means that the semivariogram starts at a value greater than zero. The nugget effect can be included in any of the semivariogram models by adding a constant. 156

Mapping Properties and Uncertainties

50-

jz

40-

0

o

S 30(0

1(D

CO

Sill = 26

o

o on ^"

O

O

/o

10-

yo Range = 6

0-c

1

' 2

' 3

' 4

' 5

' €;

' 7 ' 8 ' 9

'10 '

Lag Figure 7.13. Spherical semivariogram model with a sill of 26 and a range of 6 fitted to experimental semivariogram of Figure 7.11. Figure 7.13 shows a spherical model fitted to the experimental semivariogram of Figure 7.11. The model has a range of 6 lags, equal to a distance of almost 4000 ft. The sill has been set to 26, which is slightly greater than the variance of the 20 observations along the line of wells. Although calculation of the experimental semivariogram ideally involves use of observations equally spaced along a straight traverse across the surface, conditions are seldom ideal. In some circumstances it may be possible to devise a more or less straight path from well to well across a densely drilled area, constructing a traverse along which the observations are approximately equally spaced, and so compute a semivariogram in the classical manner. However, most people would prefer a computerized procedure for creating the semivariogram automatically, regardless of the locations and spacings between the wells. Such a procedure can be developed using the directional search capabilities of a contouring program, and is included in the RISKMAP software. An automated semivariogram calculation procedure requires that a number of approximations be made. First, since we are unlikely to find 157

Computing Risk for Oil Prospects — Chapter 7 a single, straight traverse of points in the data set, the program uses the octant search procedure to locate pairs of points that lie within the same broad, general orientation. For example, the program might find all pairs of points that are generally aligned north-by-northwest; the compass orientations of the lines connecting such pairs could vary from N to N45°E (or S to S45*^W). For the purpose of calculating a semivariogram, all of these pairs of points are considered to reflect the variation in the surface along a specific, average orientation. Similarly, the pairs of points that are found will not all be separated by the same distance; indeed, except for the effect that may be caused by a minimum allowable well spacing, in a typical set of wells we expect to see a continuous distribution of the distances between pairs of wells. This continuous range of distances is divided into discrete intervals, and all pairs of wells which fall within an interval are considered to be separated by a constant distance. It is as though we have divided the space around each well into a pattern of bins, and we pretend that any other well that falls into a bin is located at the exact center of the bin. In this way, we achieve a set of equidistant pairs of points whose squared differences can be averaged to yield an estimate of the semi variance. Of course, we've introduced inaccuracies into the estimation because the pairs of points are not really located along straight lines at equal intervals, but we hope that this uncertainty is overwhelmed by the very large number of observations that can be included if we use these approximations. Figure 7.14 shows the experimental semivariogram for thickness of the XVa Limestone as calculated by the automated procedure in RISKMAP. It is necessary to specify the radius of the bin; in this example the radius of a bin has been set so the effective distance for one lag is 1 km. The program automatically sets the radial segments at successive 45° intervals. Semivariances can be calculated for any of these preset orientations; or if there are no significant differences with orientation, they can be averaged to yield an omnidirectional semivariogram. Thickness does not show a pronounced grain, so an omnidirectional semivariogram has been calculated and is plotted in the figure. The program also will fit a semivariogram model based upon specified parameters. In Figure 7.14, a spherical model has been fitted, after some experimentation, using a sill of 21 m^ and a range of 13 km. This model appears to give a good approximation to the experimental semivariogram over a distance of about 10 km or so. Since the influence of the nearest wells is most critical, and at most locations the nearest wells will be within 10 km of the grid node being estimated, this model seems adequate.

158

Mapping Properties and Uncertainties 3025Sill = 21 o

o

o o

"~T—

15 Distance, km

—}—

20

25

Figure 7.14. Spherical semivariogram model fitted to thickness of XVa Limestone in 83 wells in the Magyarstan training area using modeling procedure in RISKMAP. The sill is 21 m^ and the range is 13 km.

Kriging The purpose of determining a model for the semivariogram is to provide information on the spatial rate of change to a procedure for estimating the form of the surface. The estimation procedure is kriging, which is a form of weighted averaging just like the procedures used in ordinary contour mapping. We can use kriging to estimate values of a surface at the nodes of a regular grid across a map, and then use contour-drawing procedures identical to those included in ordinary mapping packages to construct a contour map. The difference is that the map produced by kriging will be based on a weighting function that is tailor-made for every location in the map, and which varies with the distances and arrangements of the wells that are used in each estimate. The weights are found from the semivariances that correspond to the distances between each of the wells used in the estimate, and the distances between these wells and the point where the estimate is required. These semivariances are entered into a set of simultaneous equations, one for each well used in the estimation process. The set of equations is solved to yield the weights that are applied to each of the wells when their values are averaged to form the estimate. Actually, it is 159

Computing Risk for Oil Prospects — Chapter 7 necessary to include another equation in the set, which acts to constrain the weights so that they sum to 1.0; otherwise, the surface may be consistently biased. The kriging algorithm implemented in RISKMAP searches for the nearest well within each of the eight octants around the grid node being estimated. Therefore, a set of nine simultaneous equations must be solved and evaluated for every node in the surface grid; as you can imagine, this is computationally more demanding than ordinary contour mapping. Davis (1986) provides a simple numerical example of kriging that shows how the values from the semivariogram model are inserted into the kriging equations, how the set of equations are solved, and how the resulting weights are used to estimate the value at the location being evaluated. We will not repeat this development here, but simply show the basic equations. The simplest set of kriging equations that can be solved consists of four simultaneous equations and is based on the semivariances corresponding to the distances between three wells. The expansion of this equation set to eight wells is straightforward. ^ i 7 i i + ^2712 + ^^3713 + A = 7ip Wi^n

-h ^t^2722 -h ^/^3723 + A = 72p

w;i7i3 H- t(;2723 + ^^3733 + A = 73p

The unknown wis are the weighting coefficients to be determined, and the 7ij's are the semivariances taken from the semivariogram model. For example, 713 is the semivariance corresponding to the distance between well 1 and well 3, and 72^ is the semivariance corresponding to the distance between well 2 and the point being evaluated. The set of equations can be rewritten and set in matrix form for solution: 711

712

713

1"

712

722

723

1

7l3

1

723

1

1

733

0

1

.

'Wl' tV2 Ws

. A .

^

•71P 72p 73p

. 1

Once we've found the weights, they are used to estimate the value of the surface at locations p. We simply multiply the value Yi at each of the wells by its corresponding weight and add the products: Yp = wiYi -f ^V2Y2 + ivsYs However, we've also estimated a fourth coefficient. A, which can be used, along with the semivariances, to estimate the uncertainty in the kriged 160

Mapping Properties and Uncertainties value. This is done by weighting each of the semi variances for the distances to the point, and adding the products to A:

This value is called the error variance; conventionally, we take its square root so its units are the same as the units of measurement of the variable being mapped. Then, it is the standard deviation of the error (or simply, standard error) in our kriged estimate. For every point on the grid, we have calculated two numbers: the estimated value of the surface itself, and the standard error of that estimate. Both numerical grids can be displayed as contour maps. Figure 7.15 is a contour map of the kriged thickness of the XVa Limestone measured in the 83 wells of the Magyarstan training area data set. The mapping conventions (scale, contour interval, symbols, etc.) are the same as used in the conventionally produced contour map of thickness shown in Figure 7.2. Careful comparison of the two maps will show that near the data points the contours have the same values, but the maps convey distinctly different general appearances. Which is best for exploration purposes? Such a question can be answered only by drilling, but the map produced by kriging does have one distinct advantage that is especially significant for risk analysis. Figure 7.16 is the second map produced by kriging the thickness data; this map shows the standard error of the estimated surface. The standard error is given in the same units as the original map, or in this example, in meters of thickness. The values on the contour map give the width of one standard deviation about the estimated value of the surface. Near wells, the uncertainty about the thickness of the XVa Limestone is very low (at the exact locations of the wells we know the true thickness, so the standard error is 0.0). At a great distance from the nearest wells, such as a location in the center of the map, the uncertainty is relatively large and the standard error exceeds 3 m. If we assume that the kriging errors of estimation follow a normal distribution, we can make a probabilistic interpretation of the standard error map. Imagine that it would be possible to drill another set of 83 wells in this area, with a somewhat different configuration of locations. If we calculated a semivariogram from these new wells, it should be very similar to our original, and if we drew a new map using kriging, it should also be very similar. In fact, if we were to repeat this imaginary exercise many times, we would expect that at any location, 68 out of 100 times we would estimate a value by kriging that would lie within plus or minus one standard error of the original kriging estimate. Ninety-six times out 161

Computing Risk for Oil Prospects — Chapter 7

6490.

6485.

6480.

6475.

6470.

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.15. Isopach map produced by kriging thickness of the XVa Limestone in the Magyarstan training area. Contour interval is 2 m. Compare map to conventionally produced contour map in Figure 7.2. of 100, we would expect the estimates at a location to fall within plus or minus twice the standard error of the original kriging estimate. These expectations are derived from the properties of the normal distribution, as can be seen in Figure 7.17. A common (but arbitrary) choice for probability limits is 90%; if we set hmits of (±1.63 x standard error) around the kriging estimate, we have defined an interval that should capture the true value of the mapped surface 90% of the time. That is, the probability is only 10% that the true value at that location lies outside the specified limits. There is a 5% probability that the true value is below the lower limit, and a 5% probability that it is above the upper limit. 162

Mapping Properties and Uncertainties 6490.

6485.

6480.

6475. h

6470. r

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.16. Contour map showing magnitude in meters of one standard error in the isopach map of thickness of XVa Limestone in the Magyarstan training area. Contour interval is 0.5 m. Since the kriging operation has produced estimates of the thickness of the XVa Limestone at a regular array of locations across the map area, and estimates of the standard error at each of these locations, we can combine the two grids to create maps of the greatest likely and smallest likely thicknesses for the unit. Figure 7.18 shows a map obtained by multiplying the values in the grid of standard errors by 1.63 and then adding the result to the grid of kriged estimates of thickness. The map conventions are the same as in Figure 7.2. At the well locations, the greatest possible values of thickness are equal to the measured thicknesses in the wells, since the standard error is zero. Away from the wells, the possible thickness could be greater, because of the reduced amount of control. Figure 7.19 shows 163

Computing Risk for Oil Prospects — Chapter 7

Lower Limit

Kriging Estimate

Upper Limit

n Standard Error

- 3 - 2 - 1 0 1 2 3 I I I I I I I I I I I I I I I I I I I I I Thickness, m 20 25 30 35 40 Figure 7.17. Distribution of error around a kriged estimate of thickness of XVa Limestone in the Magyarstan training area, where the estimated thickness is 30 m and the standard error is 3 m. Upper and lower Hmits are defined as (±1.63 X standard error), forming an interval from 25 m to 35 m which contains the true thickness with 90% probability.

the corresponding lower limits on thickness, produced by subtracting (1.63 X standard error) from the kriged thickness. Perhaps the easiest way to visualize the concept of the kriging estimate and its standard error is to examine a profile across the map. Figure 7.20 shows a plot of thickness of the XVa Limestone as estimated along row 10 of the map grid. This is an east-west line at the X2-axis coordinate 6464, which passes through two dry holes in the eastern part of the area. On the profile, the kriged estimates along this row are shown by circles, each one representing a column in the grid. Near the two wells on the profile, the upper and lower limits pinch together and become equal to the kriged estimate, which in turn is equal to the measured thickness in the wells. Between the wells and elsewhere along the profile, the upper and lower limits flare out to an extent determined by the semivariance of thickness and the configuration of wells near the profile line. We expect that the true thickness hes within this envelope. The probability is only 10% that the true thickness lies outside the envelope at any location. 164

Mapping Properties and Uncertainties

6490.

6485.

6480.

6475.

6470.

6465.

6460.

6455. 1740.

1745. 1750. 1755.

1760. 1765.

1770. 1775.

Figure 7.18. Upper limit of thickness of XVa Limestone in Magyarstan training area found by adding (1.63 X standard error) to thickness estimated by kriging. The probability is 5% that the true thickness exceeds this limit.

A Complication and a Way Out Occasionally, we will encounter a geological variable which perversely refuses to yield a well-behaved semivariogram such as that shown in Figure 7.14 for the thickness of the XVa Limestone. Instead, the semivariance will increase with increasing lag, showing every indication of continuing upward forever. No sill will be apparent in the plotted semivariogram. As an example, we can analyze the data on structural elevation of the XVa Limestone in the file TRAINWEL.DAT, and we will see that its semivariogram has these unfortunate properties (Fig. 7.21). This presents serious problems for kriging, because we cannot define either a sill or a range on 165

Computing Risk for Oil Prospects — Chapter 7 6490.

6485.

6480.

6475.

6470.

6465. h

6460. h

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.19. Lower limit of thickness of XVa Limestone in Magyarstan training area found by subtracting (1.63 X standard error) from thickness estimated by kriging. The probabiHty is 5% that the true thickness is less than this limit. such a semivariogram and we will be unable to apply any of the standard semivariogram models. The behavior implies that the influence of a point extends without limits, and that the uncertainty in estimation can increase without bounds. It also creates severe mathematical problems in the solution of the kriging equations. What is causing this condition? A basic assumption in ordinary kriging is that the geologic property is stationary; that is, its average value is approximately the same everywhere. Although the mapped surface fluctuates and may be rich in features, there is no persistent tendency for it to rise or fall in value. This is not the case with structural elevation in the Magyarstan area, as we can see in 166

Mapping Properties and

~1 5

.

1 ..

10

..J, . . . . 15

.|.,,. 20

•• 1

25



Uncertainties

1

30

1

35

Grid Column Figure 7.20. Profile along row 10 of the grids used to draw Figures 7.15, 7.18, and 7.19. The line of section runs east-west at coordinate 6464 of the X2 (north-south) axis and passes through the location of two dry holes. Kriging estimates at grid nodes are indicated by open circles.

20000-

15000H CM

E o

§ 10000H

•c

E 0) CO

o o

5000 o o 0-^>-o-ii. o o

—I—

~T—

5

10

15

20

25

Distance, km Figure 7.21. Experimental semivariogram for structural elevation of the top of the XVa Limestone in the Magyarstan training area. 167

Computing Risk for Oil Prospects — Chapter 7 1000-

7504

Sill = 750

o

Range = 10; -tT— 10 15 Distance, km

o

20

25

Figure 7.22. Experimental semivariogram for second-degree trend residuals from the top of the XVa Limestone in the Magyarstan training area. A spherical semivariogram model with a sill of 750 m^ and a range of 10 km has been fitted. Figure 7.1. This tendency for a persistent change in structural elevation was emphatically shown by trend surface analysis in Figure 7.5. A geological property that exhibits a pronounced trend is said to be nonstationary, and it cannot be correctly mapped by ordinary kriging techniques. (In fact, there will be subtle distortions of the mapped surface regardless of what mapping procedure is used; positive features will tend to be shifted up-slope and negative features shifted down-slope.) If the structural surface is leveled by removing the trend (in effect, tilting the surface back to a horizontal aspect) as was done in Figure 7.6, we will find that the semivariogram of the trend residuals exhibits more acceptable behavior (Fig. 7.22). This suggests that we could solve the problem of nonstationarity by removing a trend from the data and computing the semivariogram of the stationary residuals. Kriging could be used to map the residuals, and the resulting residual map added to the trend surface map to recreate a map of the original surface. This effectively splits the mapping problem into two parts: the fitting of a deterministic nonstationary trend and the estimation 168

Mapping Properties and Uncertainties of the stationary residuals by kriging. Although such a procedure will work in many instances, it is considered inelegant by geostatistical practitioners who have devised a number of alternatives for estimating nonstationary properties and assessing the estimation errors. Perhaps the simplest of these is called "universal kriging," which combines the removal of nonstationarity and the estimation of the residuals into a single step. RISKMAP contains options that will perform universal kriging of nonstationary geologic variables. The set of kriging equations is expanded by adding two or more coefficients that must be estimated along with the kriging weights. The equations necessary to determine these additional coefficients contain powers and cross products of the geographic coordinates of the data points in the neighborhood around the location being estimated. That is, we incorporate the equation of a trend surface directly into the kriging equation and solve both sets of coefficients simultaneously. However, what we are finding is not really a trend surface in the sense that we have used the term earlier, because it does not extend throughout the map area but rather is confined to the neighborhood of the location being estimated. In fact, the trend surface is not actually computed at all; by calculating the kriging weights at the same time as the trend coefficients, the kriging weights are automatically adjusted for the presence of the trend. The universal kriging equations have the form 711

721

731

741

751

712

722

732

742

752

713

723

733

743

753

714

724

734

744

754

715

725

735

745

755

1

1

1

1

1

^11

-^12

-^13

Xi4

^15

^21

-^22

-^23

^24

^25

1 1 1 1 1 0 0 0

^11

-^21 '

'Wi'

"7ip

^12

-^22

W2

72p

-^13

-^23

ws

73p

Xi4

^24

W4

74p

-^15

-^25

W^

75p

0 0 0

0 0 0 .

A

1

^1

^Ip

L62J

-^2p

which is fearsome in appearance but is a simple extension of the ordinary kriging equations given earlier. Davis (1986) provides a relatively simple numerical example of universal kriging. Geostatisticians refer to the fitted polynomial function in universal kriging as a "drift" rather than a "trend" to distinguish it from the global function of trend surface analysis. After determining all of the coefficients, the surface is estimated by Yp = wiYi + W2Y2 + tv^Ys + W4Y4 4- wsYs The error variance (and from it, the kriging standard error) are found as before: si = wijip

+ ^t;272p + ^373p + ^474p + ^575p + A

169

Computing Risk for Oil Prospects — Chapter 7 The drift itself could be estimated and shown in the form of a contour map, but this is seldom done as it has no readily interpretable meaning. It is simply a mathematical construct useful to create local stationarity in the mapped property. Since the kriging weights and the drift coefficients are determined simultaneously, the kriging weights reflect the influence of the drift calculation. Because values of the surface at specific locations are estimated by an equation whose kriging weights reflect the drift, the effect of the drift is included in the final map grid. Figure 7.23 is a kriged map of structural elevation in the Magyarstan training area produced by universal kriging using RISKMAP. A first-order drift has been specified (RISKMAP allows choice of a first- or second-order polynomial drift). The accompanying map of the standard error of the estimated structural elevation is shown in Figure 7.24. To produce the kriged map, the user must provide the parameters of an appropriate model of the semivariogram that describes the spatial structure of the drift residuals, and this brings up an annoying problem. In order to apply universal kriging, we must provide the kriging program with several parameters, including the order of the drift and the parameters of the semivariogram of the residuals. However, the semivariogram of the residuals can be modeled only from the residuals, which must be determined by removing the drift from the variable to be mapped. But to remove the drift we must specify the neighborhood size (the range of the semivariogram) and order of the drift, which are part of the very semivariogram model that we seek! Breaking out of this circular impasse requires experimentation, first choosing an arbitrary order for the drift, arbitrarily selecting a model (linear, spherical, etc.) for the semivariogram, and setting a "convenient" neighborhood size. Using these assumed parameters, we then estimate the drift, find the residuals, and calculate the semivariogram of the residuals. If we've chosen appropriate parameters, the experimental semivariogram and the model semivariogram will coincide. If they don't, we adjust one or more of the parameters (neighborhood size, model, or order of drift) and try again. Since the drift is an arbitrary construct, we may expect that there will be several combinations that produce an acceptable fit. In general, it is desirable to have a large neighborhood (we are less apt to encounter locations where we cannot estimate the surface because there are too few points inside the neighborhood), so it is often better to specify a first-order drift and a large neighborhood rather than a second-order drift and a small neighborhood that gives equally good results. In our example, we can obtain very good initial estimates for the semivariogram parameters from the semivariogram of trend residuals shown in Figure 7.22. This suggests that a spherical model with a range of 10 km may be appropriate, and these are the parameters used in calculating the 170

Mapping Properties and Uncertainties

6490.

6485.

6480.

6475.

6470.

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.23. Structure contour map fitted by universal kriging to the top of the XVa Limestone in the Magyarstan training area. Contour interval is 25 m. Compare to conventional contour map in Figure 7.1. universal kriging map of Figure 7.23. Although the experimentation necessary to estimate the parameters of the semivariogram model is a nuisance, it is a necessary step if we wish to incorporate nonstationary geological properties into our analysis. As in other aspects of petroleum risk assessment, the time and effort spent in carefully establishing appropriate parameters and extracting relationships from data will be repaid in improved analyses and more reliable interpretations. Kriging can provide an important component in a petroleum riskassessment system. Suppose that in a certain play, thickness is believed to play a critical role in the localization of petroleum. Kriging can be used to map thickness, just as we have done using thickness data on the XVa 171

Computing Risk for Oil Prospects — Chapter 7 6490.

6485.

6480.

6475.

6470. b

6465.

6460.

6455. 1740.

1745.

1750.

1755.

1760.

1765.

1770.

1775.

Figure 7.24. Contour map showing magnitude in meters of one standard error in the structural map of XVa Limestone in the Magyarstan training area. Contour interval is 5 m. Limestone, and to predict locations where the thickness seems to be sufficient for oil accumulation. By assuming that the uncertainty follows a normal distribution, the kriging standard error allows us to express our estimate of thickness in probabilistic terms, and to specify likely upper and lower limits on its value. It provides a quantitative measure of the reliabihty of our perception of the mapped property, and directly furnishes a probability distribution in a form suitable for input into subsequent risk analyses.

172