Stochastic modelling of a fluvial reservoir: a comparative review of algorithms

Journal of Petroleum Science and Engineering 21 Ž1998. 95–121

Stochastic modelling of a fluvial reservoir: a comparative review of algorithms A.G. Journel a

a,)

, R. Gundeso b, E. Gringarten c , T. Yao

a

Geological and EnÕironmental Sciences Department, Stanford UniÕersity, Stanford, CA 94305, USA b Norsk Hydro a.s., Postboks 200, 1321 Stabekk, Vaekero, Oslo, Norway c Petroleum Engineering Department, Stanford UniÕersity, Stanford, CA 94305, USA Received 1 October 1997; revised 18 April 1998; accepted 18 April 1998

Abstract A typical characterization of a clastic reservoir starts with the modelling of its major compartments Žlayers, fault blocks. followed by the modelling of the distribution in space of its dominant facies, or facies associations. Many alternative approaches to the modelling of such facies exist. This paper proposes a comparative review of three algorithms representative of the spectrum available: Ž1. the truncated Gaussian algorithm, Ž2. a nested indicator Žbinary. simulation approach and, Ž3. an object-based algorithm. A common data set originating from a fluvial reservoir of the Statfjord formation ŽNorth Sea. is used. The impact on flow simulation of alternative facies models is evaluated through the performance of a tracer flow simulation using a newly developed streamlines-based simulator which does not require any prior upscaling. The crisp geometry obtained from object-based modelling of facies does not always result in significantly different Žor better. recoveries. q 1998 Elsevier Science B.V. All rights reserved. Keywords: reservoir; stochastic; facies; object; modelling; pixel

1. Introduction The distribution in space Žgeometry. of reservoir facies, more precisely facies associations corresponding to major flow units, is possibly the single most consequential heterogeneity for reservoir performance prediction ŽHaldorsen and Damsleth, 1990.. Another most important heterogeneity would be that induced by the geometry of faults ŽKay et al., 1993.. As long as the geometry of flow simulation blocks

)

Corresponding author.

can be adapted to major local heterogeneities ŽFarmer and Heath, 1990; Garcia et al., 1992., an accurate representation of facies geometry is critical for reservoir performance predictions. Because data is rarely enough to allow a deterministic Žsingle. representation of the subsurface heterogeneity, stochastic modelling or imaging of reservoir facies is becoming more prevalent. There exist several alternative approaches Žalgorithms. for the stochastic modelling of reservoir facies, each resulting in a different set of equiprobable stochastic images of the facies distribution in the reservoir. This paper based on a class taught at Stanford in the Spring of 1996 proposes a review of three algorithms

0920-4105r98r$19.00 q 1998 Elsevier Science B.V. All rights reserved. PII S 0 9 2 0 - 4 1 0 5 Ž 9 8 . 0 0 0 4 4 - 8

96

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

for facies geometry modelling representative of the spectrum of algorithms available; they are as follows: 1. The truncated Gaussian algorithm ŽMatheron et al., 1987.; 2. The indicator simulation algorithm ŽJournel and Alabert, 1988.; 3. An object-based simulation algorithm ŽDeutsch and Wang, 1996.. The first two are pixel-based, in that volumeŽs. occupied by each facies are builtrsimulated one pixel at a time. The last algorithm locates in space whole objects or whole sets of pixels at a time, each object having a well defined geometry Že.g., that of a curvilinear channel.. Since visual comparison of results may be deceptive—visually different geometries may not result in significantly different flow performances—the alternative facies models are compared through flow simulation Ž5-spot tracer simulation.. The flow simulator utilizes a newly-developed streamline-based algorithm to track particles along streamlines ŽBratved et al., 1992; Thiele et al., 1996.; this avoids blurring the comparison with the impact of usually poorly understood upscaling algorithms: the full resolution of each 3D facies model is here used without any upscaling. A common data set related to a major fluvial reservoir of the Statfjord formation ŽNorwegian North Sea. is used. Unfortunately, no production data were made available to this study. This paper is built around the five following sections. The Hekla data set is presented in Section 2. Elements of geostatistical exploratory data analysis ŽEDA. are presented in Section 3. Remarkably simple tools such as histograms and crossplots when applied to carefully selected variables can provide insight into the data set and influence critical modelling decisions. Right after volumetrics estimation Žan important step not covered here., comes the modelling of the distribution in space of lithofacies. In Section 4, it is shown that a combination of object-based and pixelbased techniques offers the versatility necessary to accommodate data from various types in a reservoir approaching maturity, i.e., where relatively dense conditioning by well data is available.

Once the geometry of lithofacies has been modelled, each facies is filled-in with numerical models of specific petrophysical variables Žporosity, permeability.. This step is summarized in Section 5. Some results of tracer flow simulation are reported in Section 6.

2. The Hekla data set The Hekla reservoir 1 is a section of a fluvial deposit in the Statfjord formation offshore Norway. The Hekla field consists of a system of rotated fault blocks with bed dipping 3–78 towards the southeast, see Fig. 1. The Statfjord formation within the Hekla reservoir is divided into two main reservoir zones, H 1 and H 2 , separated by an almost continuous mudstone layer Žnon pay., see the typical log plot of vertical well no. 2 given in Fig. 2. Within Hekla, the formation is an upward coarsening sequence of interbedded mudstones, siltstones and sandstones associated to braided fluvial channels with low sinuosity. The conceptual depositional fluvial model proposed to guide the numerical modelling exercise is given in Fig. 3. No vertical communication between H 1 and H 2 is observed from repeat formation testing ŽRFT data.. Lateral communication in H 1 is good, with an average net to gross ratio of 45%, an average porosity of 28%, and an average pay permeability of 400 md, see Table 1. There appears to be less lateral communication in the lower zone H 2 which has also significantly lesser reservoir qualities. The modelling presented here is limited to the upper zone H 1 , restricted to a rectangle S Ž5000 m East and 6500 m North., see Fig. 4. The data made available for this study are: Ø Well data originating from 20 wells, most deviated and drilled from two platforms, see Fig. 4. Ø Seismic data Ždepths and impedances. defined on a grid 50 m = 50 m. The surface map shown in Fig. 1 is built from the seismic depth data. 1 For confidentiality reasons, the reservoir name has been changed and its coordinates system modified Žrotation and affinity.; also petrophysical data were multiplied by constant factors. These transforms do not affect the geological structure nor the rank order statistics of the original data set.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

97

Fig. 1. Top surface of reservoir zone H 1 as interpolated from seismic data Ždepths..

The well logs include the zonation code ŽH 1 and H 2 . as defined by the geologist, the acoustic impedance log, a facies log Žsee Fig. 3 for the five facies retained., log-derived porosity and permeability, and, for the vertical wells, some core porosity and permeability Ž K V and K H . data. Fig. 2 gives a typical log plot corresponding to vertical well no. 2. Well nos. 1–4 are vertical, well no. 8 is horizontal, while all the other wells are deviated although they are near to vertical within the study area. The 3D locations of all wells are assumed accurate. Prior to the study presented hereafter, the all-important delineation of large scale reservoir compartments Žfault blocks, layers. and the estimation of corresponding volumes were done and the results are accepted here without further questioning.

3. Exploratory data analysis (EDA) Besides getting familiar with the data available, EDA typically helps into the decision of stationarity, that is, into deciding which pools of data should be retained to allow reasonable statistical mass without blurring important differences. For the Hekla reservoir, two critical decisions had to be made: Ø Should layers H 1 and H 2 be separated, Ø What are the possible facies associations regrouping some of the five original facies Žsee Fig. 3..

3.1. Stratigraphic coordinate correction Because of layer dip and uneven thickness, lateral geological continuity is calculated along ‘strata’ corresponding to a constant standardized time horizon t defined as: t s w z y z bot x w z top y z bot x g w 0,1 x

Ž 1.

where z is the vertical Ždepth. coordinate of a 3D location u s Ž x, y, z .. z top and z bot are the depths of the layer ŽH 1 or H 2 . top and bottom at horizontal location Ž x, y .. Note: z top and z bot should relate to the reconstructed layer prior to any dip, erosion or faulting. The gentle dip present over Hekla was not a concern. Similarly, the erosion of both layers tops is sufficiently limited in its areal extent to be ignored.

3.2. Facies properties The vertical proportion p k Ž t j ., for each of the five original facies k s 1, . . . ,5 and for 20 equal time intervals w t jy1 , t j x, t j s jr20, j s 1, . . . ,20, are given in Fig. 5. Zones H 1 and H 2 are separated by a continuous non-pay mudstone layer Žfacies 5.; this lack of verti-

98

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 2. Logplots of vertical well no. 2. The bottom of zone H 1 is interpreted at depth 1958 corresponding to a discontinuity in the acoustic log and the lithology log.

cal communication between H 1 and H 2 is confirmed by examination of the logs of the wells intersecting their interface and by RFT data. Rough areal facies proportions in H 1 and H 2 were obtained by kriging the facies proportions averaged over those wells with little horizontal deviation. The

high pay areas in H 1 and H 2 do not match: this implies that the geometry of the channel complexes is different in H 1 and H 2 , another reason to separate these two layers into two distinct sequences. However, the statistics of core porosity per facies type are similar across H 1 and H 2 , see Table 1. The

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

99

Fig. 3. Conceptual model for a braided fluvial channel system in the Hekla area.

corresponding facies-specific statistics, including variograms, could be averaged over H 1 and H 2 for greater statistical mass. The core K H data for facies 5 Žmudstone. appear suspect Žtoo high mean in H 1 , and too high variance in H 2 .; for mudstone the H 1 mean values were ignored. 3.3. Facies associations Facies-specific cumulative histograms Žcdf’s. of petrophysical properties were run regrouping data

from H 1 and H 2 and are shown in Fig. 6. Clearly facies 2 Žchannel sand. constitutes the best reservoir rock while facies 5 Žmudstone. is the poorest. The other sands Žconglomerate facies 1, and levies and crevasse facies 3 and 4. are intermediate in reservoir quality. If one ignores facies 1 whose proportion and spatial connectivity is insignificant, one could group facies 3 and 4 into a facies association called ‘border’ sands abutting on the channel sand facies 2. Such grouping would be consistent with the depositional model sketched in Fig. 3.

100

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Table 1 Facies proportions and statistics of core data Žnumber, mean and coefficient of variation, H 1 vs. H 2 . Facies type

Statistics

H 1 zone Prop Ž%.

1

2

3

4

5

All pay Ž1 q 2 q 3 q 4.

ns ms cv s ns ms cv s ns ms cv s ns ms cv s ns ms cv s ns ms cv s

2.2

29.7 0.19 5.8

7.5

54.8

45.3

H 2 zone

f Ž%.

KH Žmd.

12 22 0.38 235 29 1.46 31 26 0.27 52 23 0.30 62 21 0.36 330 28 0.24

12 115 1.24 234 528

3.4. Seismic impedance Only one surface seismic attribute was made available for this study, seismic impedance selected as the one best correlated with porosity. The lateral resolution of seismic impedance data offers a potential source of soft data that could help the areal mapping of facies types and related petrophysical properties. That seismic data is only 2D having no vertical resolution over the average thickness Ž25 m. of either layer H 1 or H 2 . Because of that difference in scale, the integration of seismic data into 3D mapping of petrophysical properties is difficult ŽDeutsch et al., 1996.. A 2D window was defined around each seismic impedance datum sŽ x, y ., e.g., the window could be the seismic grid cell 50 m = 50 m. All facies indicator data from well intersections of a given layer H 1 or H 2 which project horizontally onto that window were averaged together; the corresponding facies proportions p k Ž x, y ., k s 1, . . . ,5, were then plotted vs. the impedance value sŽ x, y ., figures not shown. Utilization of larger window sizes allows averaging more facies indicator data into more representative

28 894 1.95 51 167 3.41 54 168 2.46 325 488 1.81

Prop Ž%. 1.0

16.6 0.30 3.4

10.1

68.9

31.1

f Ž%.

KH Žmd.

5 26 0.35 124 27 1.63 16 21 0.39 30 22 0.35 52 14 0.39 175 25 0.33

5 223 0.83 121 517 14 163 2.55 30 72 1.83 43 39 4.69 170 401 1.86

facies proportions. All results indicate essentially a zero correlation. Many combinations of window sizes and facies associations were tried to no avail. The poor correlation of seismic impedance with facies data contrasts with the reasonable correlation between acoustic impedance and core porosity along wells and across all facies types Žlinear correlation r up to y0.57 and rank correlation r up to y0.65; figures not shown.. This indicates that there is substantial vertical facies variability within a layer ŽH 1 or H 2 ., a variability not captured by the seismic impedance data which have no vertical resolution over H 1 or H 2 . Fig. 7 gives the experimental semivariograms of seismic impedance in the two directions N308E and N608W aligned and perpendicular to the direction of channelling in layer H 1. The excellent lateral continuity displayed by the seismic data may incite to borrow parameters Žrange and anisotropy. from their variograms to model horizontal variograms of facies proportions and related petrophysical properties. We warn against such borrowing before a significant level of correlation between seismic data and well data has been established.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 4. Study zones and well locations.

Various filtering of the seismic data ŽYao and Mukerji, 1996. were performed and the results correlated to facies proportions and vertically-averaged porosity. All attempts resulted in non-significant correlation. It was thus decided to forfeit using seismic impedance for lateral resolution of facies proportions.

4. 3D modelling of facies associations The following modelling decisions were taken based on the exploratory data analysis and geological interpretation: Ø Layer H 1 and H 2 are to be modeled separately; only the results for layer H 1 are presented hereafter. Ø Three facies associations ŽFA. are retained, in lateral sequence going inward from the outer most, see Fig. 3: FA1: host sediment, i.e., facies 5 Žmudstone., FA2: ‘border’ sands, i.e., facies 3 q 4 Žlevies q crevasse splays.,

101

FA3: channel sand including the small proportion of conglomerate sand, i.e., facies 1 q 2. Prior to modelling of petrophysical properties, the distribution in the 3D space within H 1 of these three FAs must be mapped. Many alternative approaches to the modelling of categorical variables exist. The following approaches representative of the spectrum available were retained: Ø 4.1, the truncated Gaussian algorithm ŽMatheron et al., 1987., a pixel-based technique whereby a continuous 3D surface is truncated by Ž K y 1. threshold surfaces to yield K categories distributed in space. Ø 4.2, a nested indicator simulation algorithm ŽJournel and Alabert, 1988., another pixel-based technique whereby binary indicators of FAs are directly simulated in space accounting for the nested nature of their geometry. Ø 4.3, an object-based technique ŽGeorgsen and Omre, 1993; Deutsch and Wang, 1996., whereby individual channel geometries are simulated as objects. The border sands ŽFA2. are simulated independently using a pixel-based indicator algorithm and are attached to the objects ‘channel’; the main reason is that those border sands do not have the crisp geometry of channels and are less amenable to object-based simulation. Facies proportion data. All approaches considered must honor all of the following conditioning data: Ø Facies indicator data at well locations Žthese are considered hard data.. Ø Facies associations global proportions, p k for facies association k. Ø Vertical proportion curves, p k Ž t ., at time horizon t, as depicted in Fig. 5. These vertical proportions are typically obtained from averaging laterally Žalong a constant time horizon. well facies indicator data. This lateral averaging may involve declustering weight if there is preferential clustering of wells in high pay zones. Ø Lateral proportion curves, p k Ž x, y ., at horizontal location Ž x, y .. Ideally, these lateral proportions are given by seismic data, a source of information distinct from the well data. Because of the poor correlation between seismic data and vertically averaged facies proportions at well locations, this lateral information was not available for the present study.

102

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 5. Vertical facies proportions in layer H 1 Žvertical stratigraphic coordinate t ..

All previous conditioning facies proportions are corrected to be consistent one with another. Simulation grid. For the modelling of the facies associations geometries and petrophysical properties, a 3D simulation grid with mesh size 50 m east, 50 m north and 20 time intervals, including a total of 101 = 131 = 20 s 264,620 nodes was retained. This grid covers the entire area displayed in Fig. 1.

tions were considered here, with facies 1 Žconglomerate sand. separated from and nested in facies 2 Žchannel sand..

Algorithm recall. Let y l Ž u. be the lth realization of a 3D standard normal Gaussian random function Y Ž u. with covariance r Y Ž h.. Let y k Ž u. be the k th threshold surface Žsame for all realizations of Y Ž u... The facies simulated at location u is k if:

4.1. Truncated Gaussian algorithm y l Ž u . g Ž y ky 1 Ž u . , y k Ž u . x ,k s 1, . . . , K Multiple truncation of a same continuous surface, see 1D example in Fig. 8, yields nested sets well suited to represent the distribution in space of nested categories or facies. Ž K y 1. truncations would yield K categories. Since the number of truncations is not a concern with this algorithm, four facies associa-

with by construction: yo Ž u. s y`, yK Ž u. s q`, ; u. The threshold surfaces y k Ž u., k s 1, . . . , K y 1, are determined to match the facies k proportion p k Ž u. at location u s Ž x, y,t .. Note that, depending

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

103

Fig. 7. Experimental semivariograms of seismic impedance in directions N308E and N608W, along and across direction of channeling Žlayer H 1 ..

Fig. 6. Cumulative histogram of core data per facies type ŽH 1 q H 2 ., see Table 1 for statistics.

on the resolution available, p k Ž u. might be constant along a same time horizon t, i.e., p k Ž u. s p k Ž t ., or completely stationary, i.e., p k Ž u. s p k . The corresponding k th threshold is:

ŽDeutsch and Journel, 1997.. Program gtsim performs the truncations. The 3D variogram g Y Ž h. s 1 y r Y Ž h., that of Ž Y u. for separation vector h, can be synthesized from: Ø The vertical indicator variogram of the sum of pay facies Žhere 1 q 2 q 3 q 4., identical to the variogram of the complement non-pay facies Ž5, i.e., mudstone.; Ø Geologic interpretation of horizontal anisotropy and sinuosity; Ø The horizontal variogram of seismic impedance provided that a reasonable correlation has been established between seismic data and facies proportions Žnot the case here.. Reproduction of hard facies indicator data at well locations u a is obtained by conditioning the realiza-

y k Ž u . s Gy1 Ž p k Ž u . . with: Gy1 ŽP. being the inverse of the standard normal distribution function, p k Ž u. s Ý kk Xs1 s pXk Ž u. is the k th cumulative proportion. For the simulation of the continuous Gaussian field Y Ž u., many algorithms could be used. Within the public-domain software GSLIB we retained the sequential Gaussian algorithm and program sgsim

Fig. 8. Multiple truncation of a continuous field.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

104

tion  y l Ž u.,u g A4 to pseudo y-data determined from the facies data as follows: y Ž ua . g Ž y ky1 Ž ua . , y k Ž ua . x , if facies k is observed at ua , e.g., y Ž ua . s Gy1

ž

p ky 1 Ž ua . q p k Ž ua . 2

/

Results. The variogram g Y Ž h. adopted for the Gaussian random function Y Ž u. is spherical with zero nugget effect, range 7000 m in the N308E direction of channeling, 700 m in the direction N608W across channelling and four stratigraphic intervals in the vertical direction Ži.e., one fifth of the local thickness of H 1 .. The vertical proportion curves of Fig. 5 were imposed, thus defining threshold curves which do not vary laterally, y k Ž u. s y k Ž t .. The realization y l Ž u. does, though, vary in 3D. Conditioning pseudo y-data were generated at all 4396 well locations where facies data are available. Fig. 9 Žleft column. gives three stratigraphic horizons of one simulated realization of the four facies associations. The nesting of facies associations from the outward host mudstone Žfacies 5. is well reproduced, however, the FA 3 q 4 appears too scattered. The vertical facies proportion curves are correctly reproduced except for the less frequent facies 1 Žconglomerate., figures not shown. Restitution of the actual vertical coordinate by back transform of the coordinate transform Ž1. allows to map vertical sections of the simulated facies geometries showing the original dip and fault blocks, see Fig. 10. Cleaning the simulated images. The salt and pepper appearance of the simulated images of Fig. 9 Žleft column. are typical of pixel-based simulation algorithms. The noise from these images can be cleaned using various algorithms Že.g., Deutsch and Journel, 1997.. They amount to regrouping isolated pixels of any given facies while preserving its proportion and the exactitude property Žreproduction of hard data.. Fig. 9 Žright column. gives the result of such cleaning using GSLIB program trans. Note that the degree of cleaningrsmoothing can be controlled and need not be as dramatic as in Fig. 9.

Appraisal. Inspection of the results of Fig. 9 reveal the well-known drawback of pixel-based simulation algorithms in general, and of the Gaussian truncated algorithm in particular. Ø Simulated geometries are noisy, they do not present the crisp boundaries of geologist drawings. Such noise can be corrected to a certain degree. Ø Simulated geometries are rectilinear: Fig. 9 does not display any meandering of the channels Žfacies 2.. Curvilinear geometries could be obtained by pre-simulating variations in space of the anisotropy angleŽs., ŽXu, 1996.. All facies simulated by truncation of a single anisotropic surface will share that anisotropy. One would like crevasse splays Žfacies 4. and floodplains Žnot considered here. not only to be attached to channel but also to have elongation in a direction opposite of that of the channel. This is not possible with the present implementation of either programs heresim or gtsim ŽDeutsch and Journel, 1997.. Conversely, pixel-based algorithms can honor exactly any density of hard conditioning data Žwell data.. The simulated images of Fig. 9 honor exactly by construction all 4396 hard facies indicator data; this is done without discontinuities at hard data locations. Utilization of the sequence of programs sgsim and gtsim of GSLIB version 2.0 is reasonably straightforward and the cpu demand is reasonable, equivalent to any run of a sequential Gaussian simulation Žprogram sgsim alone..

4.2. Sequential indicator simulation algorithm Sequential indicator simulation ŽSIS. is a pixelbased simulation algorithm that builds a categorical image, pixel after pixel, by drawing from the categories local probability distributions ŽJournel and Alabert, 1988; Alabert and Massonat, 1990.. These probability distributions are updated to account for the categories already simulated at neighboring pixels. Characteristic of pixel-based algorithms, the resulting geometries Žfor each category. tend to be noisy; conversely, honoring exactly dense conditioning data is not a problem. Another advantage of indicator techniques is their flexibility in accommo-

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

105

Fig. 9. Three stratigraphic horizontal levels of a 3D simulation of four facies associations by multiple truncation of a simulated 3D continuous surface. Original simulation Žleft column. and cleaned images Žright column..

106

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 10. A vertical section stricking N608W of the facies geometry 3D simulation after actual vertical coordinate restitution. Vertical exaggeration is 5r1. ŽNote the continuity conformable to top and bottom surfaces and the three fault blocks..

dating different indicator variograms for different facies indicators and, more importantly, in accommodating diverse types of soft data ŽZhu and Journel, 1993; Deutsch and Journel, 1997.. Because of the nesting of the three facies associations retained Žmudstone, border sands, channel sands., a sequence of two nested binary indicator simulations was considered instead of a single multiple Žthree categories. indicator simulation. Nested indicator simulation. The following sequence was adopted: Ø Step 1: Simulate a binary indicator separating the channel sands FA Žchannelq conglomerate. from the rest. Ø Step 2: Within the domain Žpixels. simulated as non-channel sands, simulate a second binary indicator separating the border sands FA Žlevies q crevasses. from the mudstone. Proper care is given to ensure that the simulated border sands be ‘attached’ to the previously simulated channel sands, see the conceptual model of Fig. 3. Step 1: Channel vs. non-channel Let I1Ž u. be the first binary indicator defined as: I1 Ž u . s 1 if u g channel sands, s 0 if not A sequence of indicator krigings provides at each node the probability that I1Ž u. s 1, conditional to neighboring original well indicator data and previ-

ously simulated I1-values. Let p 1 )Ž u. g w0,1x be that conditional probability. A random number Õ uniformly distributed in w0,1x is drawn, then: if Õ F p 1Ž u., pixel u is assigned to channel sands: i 1Ž u. s 1, if Õ ) p 1Ž u., pixel u is assigned to non-channel sands: i 1Ž u. s 0. The simulated indicator value i 1Ž u. is then used as datum for the determination of the conditioning probability at nodes uX neighbors of u. The indicator kriging of p 1Ž u. requires the 3D variogram of the indicator variable I1Ž u.. Vertical well data provides a reasonably clear experimental vertical variogram with a range of about four stratigraphic time units Žout of a total of 20., see Fig. 11a. There is not enough lateral well data to expect any interpretable horizontal variogram Ževen including the short span of horizontal well no. 8.. The poor correlation between seismic impedance and facies proportions prevents borrowing range parameters from the seismic impedance variograms ŽFig. 7.. This is when geological interpretation becomes essential. It was considered that the average channel complex width is about 500 m leading to retaining a 500-m range in the direction N608W. Along the N308E direction of channeling, the range was assumed to be 5000 m, for a horizontal anisotropy ratio of 1:10. An anisotropic spherical variogram with zero nugget effect was used; the sill plays no role and was set to 1. The simulation was done in the stratigraphic coordinate system Žwith 20 vertical time intervals. and the results were reinstated into the original coordinate system.

Vertical facies proportions. The vertical channel sands proportions p 1Ž t ., t s 1, . . . ,20, as deduced from well data Žsee Fig. 5., are accounted for by considering a simple kriging system with a vertically varying mean. The estimated local probability for I1Ž u. to be 1 Žchannel sands. at location u s Ž x, y,t . is written, ŽDeutsch and Journel, 1997, p. 64.: n 1Ž u .

p1 ) Ž u . s p1Ž t . q

Ý as1

la i 1 Ž ua . y p 1 Ž ta .

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

where: i 1Ž ua . is the well indicator data at location ua s Ž x a , ya ,ta ., a s 1, . . . ,n1Ž u., n1Ž u. is the number of facies indicator data found in the neighborhood of u. p 1Ž t . and p 1Ž ta . are the prior vertical

107

proportions of border sands at time horizons t and ta , la is the simple kriging weight. Note that at any datum location u s ua 0 , the kriging system yields la 0 s 1, la s 0, ;a / a 0 ,

Fig. 11. Vertical facies indicator variograms in layer H 1 .

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

108

thus: p 1 )Ž u. s i 1Ž ua 0 . equal to 0 or 1, hence the simulated value i 1Žu. drawn at u necessarily identifies the indicator datum value i 1Ž ua 0 .. Step 2: Border sands within non-channel sands This second step is limited to Žnested into. those nodes u that have been simulated as being non-channel sands, i.e., i 1Ž u. s 0. Consider a second binary indicator defined as: I2 Ž u . s 1 if u g border sands Ž levees q crevasse. , s 0 if not The simulation algorithm for I2 Ž u. is in all points similar to that presented in Step 1 for I1Ž u. with the notable addition of an attraction parameter to ‘attach’ the simulated border sands, i 2 Ž u. s 1, to previously simulated channel sands. Attracting border sands to channel sands. The local probability for I2 Ž u. to be 1 Žborder sands. is written: n 2Ž u .

p2 ) Ž u . s p2 Ž t . aŽ u . q

Ý

na i 2 Ž ua .

as1

yp 2 Ž ta . where: i 2 Ž ua . is the well indicator data at location ua s Ž x a , ya ,ta ., a s 1, . . . ,n 2 Ž u., n 2 Ž u. is the number of facies indicator data found in the neighborhood of u. p 2 Ž t . and p 2 Ž ta . are the prior vertical proportions of border sands relative to non-channel sands, at time horizons t and ta . The only difference with the expression used in Step 1 is the attraction parameter aŽ u. which increases the prior probability p 2 Ž t . if location u is close to previously simulated channel sands nodes and, conversely, decreases that prior probability if u is far away from channel sands; that decrease is such that the overall proportion of border sands p 2 Ž t . at time horizon t be preserved. Variograms. The experimental vertical variogram of indicator of border sands Žwithin non-channel sands. is shown in Fig. 11b. Note the range of about two time intervals Žout of a total 20 time horizons., significantly shorter than the range 4 read on Fig. 11a for channel sands. One does expect levees and crevasses to be thinner than channel beds.

Again the horizontal variogram for I2 Ž u. had to be synthesized from geological interpretation. Border sands, in particular crevasse splays are expected to have a rounder Žmore isotropic. geometry than channels. The previously adopted anisotropy ratio of 10 to 1 for channel sands has been reduced to 2 to 1 for border sands, with 250-m range Žhalf the width of channel beds. in the N608W direction and 500 m in the N308E direction of channeling. Again, an anisotropic spherical variogram model with zero nugget effect was retained. Results. Fig. 12 gives five stratigraphic horizontal levels of the simulated border sands within the volume simulated as non-channel sands. Note the more rounded geometry of the simulated border sands and how they are attached Žin 3D. to channel sands. Fig. 13 indicates a correct reproduction of the well-derived vertical facies proportions. Ideally, results from a whole series of realizations should be shown with their cloud containing the target proportions. Appraisal. The results of Fig. 12 display again the drawback of pixel-based simulation algorithms: lack of crisp facies boundaries. Similarly to what has been done in Fig. 9, one could have ‘cleaned’ the realization of Fig. 12. The result of the nested indicator simulation approach are deemed better than that of the truncated Gaussian algorithm, compare Figs. 9 and 12, because of the added control on the geometry of the border sands given by the additional variogram of I2 Ž u.. The channel beds simulated in Fig. 12 do not undulate, the reason being the long horizontal range Ž5000 m. adopted in direction N308E, a range almost as large as the corresponding dimension of the field being simulated. Curvilinear channel sands geometry could have been simulated by considering a shorter range and spatially variable directions of anisotropy ŽAlabert and Massonat, 1990; Xu, 1996.. On some horizontal levels of Fig. 12 the border sands appear isolated, see for example the upper level t s 20. Actually these border sands are attached to channel sands previously generated at either lower or higher levels. One could control the vertical Žalong t . distribution of border sands by defining an attraction parameter that depends only on

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

the proximity of channel sands on the same level, i.e., defining a ‘horizontal’ rather than a 3D attraction parameter.

109

Running twice the GSLIB bf sisim program Žfor the two nested indicator simulations. is straightforward. Because each run of the sisim program in-

Fig. 12. Five stratigraphic horizontal levels of a 3D simulation of three facies associations in H 1 Žnested sisim algorithm..

110

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

volves only one indicator, i.e., one kriging system per simulated node, the total cpu cost is about the same as for program sgsim q gtsim Žtruncated Gaussian algorithm.. 4.3. Object-based simulation algorithm Object-based simulation algorithms amount to build in space whole objects with their specific

geometries, one at a time. Objects are added until the target number or volume proportion Že.g., net-to-gross ratio. is approximately matched. Matching other constraints such as vertical or lateral prior proportions and, more critically, local facies data at well locations, is done by adapting object densities to the local proportions and, through a trial-and-error process, modifying object sizes, moving them around, addingrremoving objects. This trial-and-error pro-

Fig. 13. Vertical facies proportions reproduction in H 1 Žnested sisim algorithm..

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

cess is typically controlled by an objective function measuring how close each constraint is approximated. The clear advantage of object-based simulation algorithms is the crisp reproduction of the target geometries. Because objects are built one at a time, a system of coordinates adapted to the specific geometry of each object can be established which is then used to simulate the petrophysical properties of that object. For example, that system of coordinates would match the sinuosity of an undulating channel Žthe object. resulting in similar sinuous petrophysical continuity. The drawback of object-based algorithms is the difficulty to honor a combination of prior proportion curves and dense well facies indicator data; this may require many trial-and-error iterations and related high cpu costs. Typically, it is difficult to honor exactly all 3D well data when the density of wells increase Žmature reservoir.. The literature on object-based simulation algorithms is extensive and their applications to simulation of reservoir facies are dominated by the Norwegian school, e.g., Bridge and Leeder, 1979; Clementsen et al., 1990; Georgsen and Omre, 1993. This section presents results obtained from the public-domain code fluvsim. The fluzsim code (Deutsch and Wang, 1996). The elongated geometry of fluvial channel complexes Žbelts. and individual channels, see Fig. 3, allows their simulation Žconstruction. as a series of very continuous orthogonal deviations from a straight center line, see Fig. 14a. Such deviations could be realizations of a continuous 1D random function; amplitude and period of the channel sinuosity are controlled by the sill and range of that random function covariance. Simulation of the fluvial deposition sequence is done one step at a time, progressing from the large scale heterogeneity Žlayers. down to the smallest scale Žpetrophysical variability within each specific channel. through a succession of reversible coordinates transform, see Fig. 14b. An iterative sequence of operations Žaddrremove, translaterrotate. controlled by an objective function allows approximation of prior constraints. The constraints allowed by fluvsim are presently: 1. Facies proportionsŽglobal, vertical andror lateral., 2. Facies indicator data Žat well locations..

111

The 1996 fluvsim code allows simulation of only one facies Žchannel sands. within its host rock Žhere mudstone.. Simulation of border sands Žlevees and crevasse splays. abutting on the channels had to be done using a sequential single indicator simulation algorithm as described in Step 2 of Section 4.2.

Results. Fig. 15 gives five stratigraphic horizontal levels of fluvsim-generated channels. Six vertical sections of the same 3D realizations are shown in Fig. 16. Conditioning includes: Ži. the global proportion: 31.9% of channel sands in H 1 , see Table 1; Žii. the corresponding vertical proportions, p 1Ž t ., t s 1, . . . 20, see Fig. 5; Žiii. all facies indicator data available from wells: there were a total of 4542 such indicator data i 1Ž ua ., mostly clustered in the center of the volume simulated. For the purpose of checking data exactitude, these well data were relocated to the nearest simulation grid Ž101 = 131 = 20.. Such relocation is not a requirement of the program. As seen from Fig. 15, no channel has been simulated on the top level of H 1 Ž t s 20., a result consistent with the input proportion of channel sands at that level. Reproduction of the vertical proportion curve, figure not shown, is poorer than that provided by the pixel-based algorithms but deemed reasonable given the expected fluctuations from one realization to another. The six EW Ž y s constant. vertical cross-sections shown in Fig. 16 also display the conditioning well facies data available on these sections: most of these well data are matched by the simulated facies without local discontinuity. Globally, there is a 85% match of the 4542 conditioning well facies data. Most of the mismatches involve only one simulation grid cell displacement. Fig. 17 compares the experimental channel sands indicator variogram Ždeduced from vertical wells. to the vertical variogram of simulated channel sands indicator values. Again the match is reasonable. Within the volume previously simulated as nonchannel sands Ž i 1Ž u. s 0., the pixel-based indicator simulation of border sands Ž i 2 Ž u. s 1. described in Step 2 of Section 4.2 was implemented. Introduction of an attraction factor allows generating the border sands next to the nodes already simulated as channel sands.

112

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 18 gives the five stratigraphic horizontal levels selected for Fig. 15. Note how the simulated border sands abut on the channels simulated by fluvsim. Reproduction of the three facies associations vertical proportions is reasonable Žfigure not shown.. Appraisal. The simulation shown on Fig. 18 results from a hybrid object-based and pixel-based simula-

tion algorithm: Ži. an object-based algorithm Žfluvsim. is better suited to simulate the elongated and sinuous geometry of channels. Note though that this algorithm is actually build from pixel-based simulation of channel center line deviations from a straight line Žsee Fig. 14a.; Žii. a pixel-based algorithm Žsisim. suffices to simulate the uncharacteristic shapes of border sands, while attaching them to the previously simulated channels.

Fig. 14. Fluvial channel simulation Žfluvsim code.. Ža. The channel sinuous center line is simulated as a series of deviations from a straight line. Channel cross-section geometry varies along its center line. Žb. Succession of coordinate transformations from layers to individual channels.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 15. Five stratigraphic horizontal levels of a 3D realization of channel sands Žlayer H 1 . as generated by fluvsim.

113

114

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 16. Six EW vertical sections of the 3D realization shown in Figure 15. The conditioning well data are shown with different colors for comparison.

The result of that hybrid simulation algorithm are deemed better than those obtained from the two strictly pixel-based algorithms, compare Fig. 18 Žfluvsimq sisim. to Fig. 9 Žgtsim. and Fig. 12 Žnested sisim.. This is, however, a visual judgement based on the crisp and sinuous geometry of the simulated channels. It is unclear that flow path in the subsurface follows such crisp boundaries and whether such visual differences make a difference in terms of oil production, see later comparison in Section 6. Also, for a more relevant comparison, the simulation of Fig. 12 Žnested sisim. should have been run with a smaller indicator variogram range in the direction N608W Žperhaps 200 m instead of 500 m. to better match the corresponding channel width of Fig. 18 Žfluvsim.. Letting program fluvsim run long enough to reach an acceptable match of the various constraints Žonly 15% mismatch of well data. took significantly Žfactor 3. more cpu than the corresponding two runs of

program sisim with two nested indicator variables. However, in both cases and for the size of the present simulation grid cpu time is not a concern.

Fig. 17. Comparison of vertical variograms, experimental from well data, and from simulated channel sands indicator values.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Fig. 18. Five stratigraphic horizontal levels of a 3D simulation of three facies associations in H 1 Žfluvsimq sisim algorithm..

115

116

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

5. Modelling petrophysical properties Within the geometric templates generated by each simulated realization of the three facies associations retained Žchannel sands, border sands and mudstone., simulation of the petrophysical properties specific to each facies associations can now be performed. Three properties were considered: porosity, horizontal and vertical permeability. The spatial distributions of these properties were modeled then simulated for channel sands and border sands, one facies association independently of the other. For mudstone

Žnon-pay facies., the properties were set to constant values taken from the corresponding core data statistics, see Table 1. Note the relatively high values for both porosity and permeability in mudstone. 5.1. Co-located cosimulation Joint simulation in space of all three variables, f , K H and K V , requires prior modelling of the joint cross-covariance matrix which includes at least six covariance and cross-covariance functions. Given the shortage of lateral Žhorizontal. information, such

Fig. 19. One realization of porosity Žlevel t s 10. conditional to a previously simulated geometric template Žthree facies associations..

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

modelling was not attempted, and the much lighter co-located cosimulation alternative was considered, ŽAlmeida and Journel, 1994.. Because it is typically better sampled and its spatial variability is both smoother and easier to model, the variable porosity is simulated first, conditional to porosity data specific to the facies association being considered. Next, horizontal permeability K H is simulated conditional to its own neighborhood data and to the sole co-located previously simulated porosity value Ždifferent from one realization to another.. Last, vertical permeability K V is simulated conditional to its own neighborhood data and to the sole co-located previously simulated K H value. This sequence of conditioning allows reproduction of the co-located cross-correlation Žat h s 0. between the three petrophysical variables. For the simulation of porosity Ža reasonably wellbehaved variable when limited to an homogeneous facies association., a multivariate Gaussian model after normal score transform was considered. More precisely the sequential Gaussian algorithm, program sgsim of GSLIB ŽDeutsch and Journel, 1997., was used. For the simulation of permeabilities, K H and K V , again, a multi-Gaussian model was considered.

5.2. The cookie-cutter approach Because pixel-based simulation algorithms Ževen in 3D as in this study. are so fast, it was deemed unnecessary to limit the simulation grid for each facies association to its previously simulated geometric template. Note that it is a trivial matter to limit simulations to a previous template. Instead, the simulation of porosity for channel sands was performed independently of that for border sands, each over the entire 3D simulation grid Ž101 = 131 = 20.; then simulated porosity specific to each facies association was retrieved Žcookie-cut. from its own simulated template and re-assembled into the final simulated porosity realization; recall that grid nodes falling within the simulated mudstone facies were assigned a constant value. Fig. 19a–c display the simulated porosity specific to each facies association simulated template; these independent realizations are re-assembled in Fig. 19d.

117

Since the petrophysical properties are here simulated independently of their template geometries, the spatial continuity of the resulting simulated values need not conform to the local geometry of the corresponding facies association. For example, the lateral continuity of simulated porosity will not meander conformably to the simulated channel sands geometry, see Fig. 19a. Object-based algorithms have the distinct advantage of yielding a system of curvilinear coordinates specific to each channel generated; conformable simulation of petrophysical properties can be achieved using that specific system of coordinates independently for each simulated channel. At the later time of modelling the second Žlower. layer H 2 , an improved version of fluvsim was available allowing channel-conformable simulation of petrophysical properties.

6. Flow simulation The final goal of detailed reservoir numerical modelling is predictionrmanagement of future performance through, e.g., flow simulation ŽAziz and Settari, 1979.. This final section is restricted to the study of the sensitivity of tracer flow responses to diverse geometries of the facies associations, for a common realization of the petrophysical variables. We also took advantage of this dense simulated data set to check the performance of the newly developed, superfast, streamline flow simulation algorithm ŽThiele et al., 1996.. 6.1. Tracer flow simulation A geostatistical numerical model is provided by any one realization of the sequence of simulations described in the previous Sections 4 and 20. The 3D simulation grid comprises 101 = 131 = 20 cells, each cell being informed with a facies type Žchannel sands, border sands, mudstone. and three petrophysical values Žporosity, horizontal and vertical permeability.. The flow simulation pattern consists of an injector located at the center of the numerical model and four producers at its four corners, a typical 5-spot pattern.

118

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Tracer is injected at the constant pressure of 4000 psi. Production at each of the four wells is set at the constant rate of 5000 m3rday. An analytical mapping of tracer flow solution along the streamlines is used. The response curves retained are the proportion of injected tracer in the total produced fluid Žreferred to as tracer-cut. and the amount Žin pore volumes. of fluid initially in place recovered Žreferred to as recovery efficiency., as functions of pore volume injected Ždimensionless time t D . for the whole field and at each of the four producing wells, see Fig. 21. Recall that these are dimensionless responses. 6.2. SensitiÕity of flow performance to channel geometry A limited study of sensitivity of flow performance to uncertainty in the geostatistical numerical model was conducted. The largest uncertainty was deemed to be associated with the geometry of the facies associations. We considered two simulated modelsrrealizations of that geometry, one generated with the pixel-based algorithm Žnested sisim, see

Fig. 21. Proportion of injected tracer in total recovered fluid at each of the four producing wells. Mudstone permeability is 39 md.

Fig. 20. Recovery of initial fluid in place and proportion of injected tracer in total recovered fluid for all four producing wells. Mudstone permeability is 39 md.

Fig. 22. Recovery of initial fluid in place and proportion of injected tracer in total recovered fluid for all four producing wells. Mudstone permeability is 0.0001 md.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

Section 4.2; Fig. 12., and the second generated with the object-based algorithm Žfluvsimq sisim, see Section 4.3; Fig. 18.. Both geometric models are conditioned to the same well facies data. To freeze the contribution of petrophysical variability, we overlayed onto these two different geometric templates the same realizations Žone per facies association. of simulated petrophysical values Ž f , K H . Žwe consider here isotropic permeability, i.e., K H s K V .. Then the fast, streamline-based, flow simulator was run with full resolution Žno upscaling. on the two numerical models. The results are shown on Fig. 20. The dimensionless responses are very similar. One possible reason is that FA5 Žmudstone. has a relatively high permeability set to 39 md. Fig. 21 shows the tracer-cut responses at each of the four wells. Tracer breakthrough occurs earlier at wells 1 and 3 which are aligned with the channel orientation. To investigate further the influence of the template geometry, a permeability of 0.0001 md was assigned to the mudstone pixels. The mudstone porosity was kept at 14% to consider the same total pore volume as the previous case. The total field

Fig. 23. Proportion of injected tracer in total recovered fluid at each of the four producing wells. Mudstone permeability is 0.0001 md.

119

responses are shown on Fig. 22. The difference between the two templates is more visible, and a later tracer breakthrough can be observed for the template generated by the object-based approach. Fig. 23 shows the tracer-cut responses at each of the four wells. The differences between the responses of each well is very clear for the template generated by the nested sisim approach. Warning: In the absence of a reference ‘truth’, no conclusion should be drawn from the limited results shown in Figs. 20–23. Also these figures give results corresponding to only one realization of each approach: in presence of sparse well data, there are large differences between realizations of a same approach, hence one could possibly generate a pair of fluvsim and sisim templates which would yield similar or very different flow performances. Time limitations did not allow developing further this sensitivity study.

7. Conclusions One main conclusion of this review study was that reservoir modelling cannot be made entirely automatic. The specificity of each reservoir geology and of each data set available requires great flexibility of the modelling tools available. There is no unique or ‘best’ approach that would be valid for all reservoirs of all types. The approach chosen to address a particular modelling problem should not be Žover. constrained by tool availability, rather, the toolbox Žsoftware. available should offer the flexibility to custom-design the algorithm deemed best suited. This calls for the reservoir modeler to be conversant with concepts and algorithm development. No matter the experience of the modeler, that experience may not be tapped to its full extent if constrained by the limitations of a rigid software. The application of geostatistics, in particular, and the art Žin the best sense of the term. of reservoir modelling, in general, call for a series of extremely consequential yet somewhat subjective early modelling decisions. Examples are decisions about stationarity pools, whether of geological facies or petrophysical populations, decisions about which data

120

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121

should be considered hard and serve as a reference against which softer data are calibrated, final decisions about which algorithms to retain when Žas often the case. the poor resolution of flow simulations cannot arbitrate. Experience is the final decision-maker. For the modelling of the Hekla reservoir it was found that a combination of object-based techniques and pixel-based algorithms offers the necessary flexibility in providing the clean channel shapes together with more random geometry of crevasse splays. Pixel-based algorithms are then used to fill the simulated geometries with petrophysical properties.

FA G h i I k K KH KV pk r s t u x y y Ž u. Y z z bot

pk g l, n f r

True vertical depth of top of H 1 layer Cumulative proportion of facies k Semi-variogram Kriging weight Porosity Correlation coefficient

Acknowledgements We thank the management of Saga Petroleum and Norsk Hydro for releasing the Hekla data set, providing the opportunity of invaluable training of several classes of graduate students at Stanford.

8. Nomenclature

aŽ u.

z top

Attraction parameter at location u Facies association Standard normal cumulative distribution function Separation distance vector Realization of binary indicator random function Binary indicator random function Facies type k Total number of facies Horizontal permeability Vertical permeability Proportion of facies of type k Rank correlation coefficient Seismic impedance Stratigraphic time horizon Location u s Ž x, y,t . Horizontal coordinate along x-direction Horizontal coordinate along y-direction Realization of random function Y at location u Standard normal random function True vertical depth True vertical depth of bottom of H 1 layer

References Alabert, F., Massonat, G., 1990. Heterogeneity in a Complex Turbiditic Reservoir: Stochastic Modelling of Facies and Petrophysical Variability, presented at 65th SPE Annual Technical Conference and Exhibition, pp. 775–790. Almeida, A., Journel, A.G., 1994. Joint simulation of multiple variables with a Markov-type coregionalization model. J. Math. Geol. 26 Ž5., 565–588. Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Elsevier, London, 476 p. Bratved, F., Bratved, K., Buchholz, C.F., Holden, L., Riesbro, N.H., 1992. A new front-tracking method for reservoir simulation. SPE RE 7 Ž2., 107–116. Bridge, J., Leeder, M., 1979. A simulation model of alluvial stratigraphy. Sedimentology 26 Ž5., 617–644. Clementsen, T. et al., 1990. A computer program for evaluation of fluvial reservoirs, In: Buller, A. et al. ŽEds.., North Sea Oil and Gas Reservoirs II. Graham & Trotman, London, pp. 373–385. Deutsch, C.V., Journel, A.G., 1997. GSLIB: Geostatistical Software Library and User’s Guide, 2nd edn. Oxford Univ. Press, New York, 368 p. Deutsch, C.V., Wang, L., 1996. Hierarchical object-based modelling of fluvial reservoirs. J. Math. Geol. 28 Ž7., 857–880. Deutsch, C.V., Srinivasan, S., Mo, Y., 1996. Geostatistical reservoir modelling accounting for precision and scale of seismic data. SPE paper 36497, Presented at the 1996 SPE Annual Technical Conference and Exhibition, Denver, CO, pp. 9–19. Farmer, Heath, 1990. Garcia, M.H., Journel, A.G., Aziz, K., 1992. Automatic grid generation for modelling reservoir heterogeneities. SPE RE 7 Ž2., 278–284.

A.G. Journel et al.r Journal of Petroleum Science and Engineering 21 (1998) 95–121 Georgsen, F., Omre, H., 1993. Combining fiber processes and Gaussian random functions for modelling fluvial reservoirs. In: Soares ŽEd.., Proceedings of the 4th International Geostatistics Congress. Kluwer, Dordrecht, 1993, pp. 425–440. Haldorsen, H.H., Damsleth, R., 1990. Stochastic modelling. J. Petrol. Technol. 42 Ž4., 404–412. Journel, A.G., Alabert, F., 1988. Focusing on spatial connectivity of extreme-valued attributes: stochastic indicator models of reservoir heterogeneities, paper SPE 18234. Kay, E.A., Bodnar, D.A., Cazier, E.C., 1993. Faulting: a major control on fluid flow and production performance, Proceedings of the 1993 Western Regional SPE Meeting, Prudhoe Bay field, AK, pp. 308–320.

121

Thiele, M., Batycky, R., Blunt, M., Orr, F., 1996. Simulating flowing heterogeneous media using streamtubes and streamlines. SPE RE 10 Ž1., 5–12. Xu, W., 1996. Conditional curvilinear stochastic simulation using pixel-based algorithms. J. Math. Geol. 28 Ž7., 937–950. Yao, T., Mukerji, T., 1996. Application of factorial kriging to improve seismic data integration. In: Baafi, Schofield ŽEds.., Geostatistics Wollongong 96, Vol. 1. Kluwer, pp. 350–361. Zhu, H., Journel, A.G., 1993. Formatting and integrating soft data: stochastic imaging via the Markov–Bayes algorithm. Soares ŽEd.., Proceedings of the 4th International Geostatistics Congress, Kluwer, Dordrecht, pp. 1–12.