ALLUVSIM: A program for event-based stochastic modeling of fluvial depositional systems

ARTICLE IN PRESS Computers & Geosciences 35 (2009) 1671–1685

Contents lists available at ScienceDirect

Computers & Geosciences journal homepage: www.elsevier.com/locate/cageo

ALLUVSIM: A program for event-based stochastic modeling of fluvial depositional systems$ M.J. Pyrcz a,, J.B. Boisvert b, C.V. Deutsch b a b

Quantitative Stratigraphy, Chevron Energy Technology Company, 1500 Louisiana St., Houston, TX 77002, USA Department of Civil & Environmental Engineering, University of Alberta, 3-133 Markin/CNRL Natural Resources Engineering Facility, Alberta, Canada T6G 2W2

a r t i c l e in fo

abstract

Article history: Received 7 February 2008 Received in revised form 27 July 2008 Accepted 26 September 2008

This paper presents an algorithm for the construction of event-based fluvial models. The event-based approach may be applied to construct stochastic pseudo-process-based fluvial models for a variety of fluvial styles with conditioning to sparse well data (1–5 wells) and areal and vertical trends. The initial models are generated by placing large-scale features, such as channels and crevasse splays, into the model as geometric objects. These large-scale features are controlled by geometric input parameters provided by the user and are placed into the model to roughly honor well data through a rejection and updating method. Yet, some model to well data mismatch may still occur due to inconsistency in the size and positioning of complicated features relative to the well data. An image processing algorithm is used to post-process realizations to exactly honor all well data. The final, cell-based models, have no data mismatch and contain geologically realistic fluvial features that would be difficult to obtain with other pixel-based methods and precise conditioning that is difficult to obtain with objectbased methods. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Reservoir characterization Geostatistics Facies modeling Object-based modeling

1. Introduction Interest in North Sea fluvial reservoirs led to the development of object-based models for fluvial facies and geometries (see Deutsch and Wang, 1996 and Patterson et al., 2002 for review of previous methods). For these models conditioning is often problematic. Difficulties in conditioning spurred research in direct object-based modeling. Viseur et al. (1998) and Shmaryan and Deutsch (1999) developed methods to simulate fluvial objectbased models that directly honor well data. These algorithms segment the well data into unique channel and non-channel facies and then fit channels through the segments. The channel centerline is parameterized as a random function of departure along a vector and the geometry is based on a set of sections fit along the center line. These techniques are well suited to paleo-valley (PV) reservoir types. The PV reservoir type geologic model is based on ribbon sand bodies from typically low net-to-gross systems with primary reservoir quality encountered in sinuous to straight channels and secondary reservoir rock based on levees (LVs) and crevasse splays (CSs) embedded in overbank fines (Collinson, 1996; Galloway and Hobday, 1996 and Miall, 1996).

$

Code available from server at http://www.iamg.org/CGEditor/index.htm.

 Corresponding author. Tel.: +1832 854 7016; fax: +1832 854 7072.

E-mail address: [email protected] (M.J. Pyrcz). 0098-3004/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.cageo.2008.09.012

More complicated channel belt (CB) fluvial reservoir types are common. Important examples include the McMurray Formation (Mossop and Flach, 1983 and Thomas et al., 1987) and Daqing Oil Field, China (Jin et al. 1985 and Thomas et al., 1987). These reservoirs include complicated architectural element configurations developed during meander migration punctuated by avulsion events. The application of the bank retreat model for realistic channel meander migration has been proposed by Howard (1992), applied by Sun et al. (1996) and Lopez et al. (2001) to construct realistic models of CB-type fluvial reservoirs; however, these methods lack flexibility in conditioning. The event-based approach (Pyrcz, 2004; Pyrcz and Deutsch, 2005 and Pyrcz and Strebelle, 2006) possesses: (1) improved flexibility to reproduce a variety of fluvial reservoir styles with realistic channel morphologies, avulsion and meander migration and (2) the ability to be efficiently updated to honor well data and areal reservoir quality trends. A robust conditioning method is introduced for settings with 1–5 near vertical wells and seismic information that identifies areal trends in net-to-gross. The associated algorithm, ALLUVSIM, includes a set of research quality building blocks for the application of event-based simulation to fluvial depositional settings with robust conditioning. This code may be applied to variety fluvial settings or adapted to other depositional settings, such as confined and weakly confined deepwater reservoirs (Pyrcz and Strebelle, 2008). This work was inspired by the developments of Sun et al. (1996) and Lopez et al. (2001), but it was conducted independent

ARTICLE IN PRESS 1672

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

of Cojan et al. (2004). The reader is referred to this work for additional insights into the construction of geostatistical pseudoprocess-based fluvial models.

2. Event-based stochastic fluvial model Event-based models are forward, process mimicking geostatistical models (Pyrcz, 2004; Pyrcz and Deutsch, 2005 and Pyrcz and Strebelle, 2006). The basic building block of event-based models is the centerline. A centerline represents the central axis of a flow event and backbone for architectural elements (Wietzerbin and Mallet, 1993). This concept is general and may represent confined or unconfined, fluvial, turbidity or debris flows. Genetically related centerlines may be grouped into centerline complexes. Centerline complexes are interrelated by process; for example, a centerline complex may represent a braided stream related by avulsion or a point bar related by meander migration. Fluvial architectural elements such as crevasse splays, channels and levees are attached to centerlines and element interrelationships are characterized by centerline complexes. This is a logical technique for constructing fluvial models, since all architectural elements are related to ‘‘flow events’’. A conceptual geometric model of the associated architectural elements is shown in Fig. 1.

3. Centerlines The application of a (random) function to represent the plan view projection of a fluvial flow event is severely limited. As a function, fs(x), may only have a single value for any value x. In graphical terms, a function may not curve back on itself. This precludes the direct use of a function (i.e., channel centerline characterized by locations x, fs(x)) to characterize high-sinuosity channel centerlines. A solution is to fit a separate cubic spline to characterize each coordinate (x and y) with respect to distance along the spline. This allows for high-sinuosity centerlines to be modeled (Pyrcz, 2004). A centerline is a table of control nodes with location and properties. Constraints must be added to prevent centerlines from crossing themselves. The advantages of this centerline technique are: (1) continuous interpolation of centerline location in Cartesian coordinates at any location along the centerline, (2) relatively few parameters required to describe complicated curvilinear paths, (3) manipulation of splines is more computationally efficient than modifying

geometries and (4) other properties such as geometric parameters and trends may be stored as continuous functions along the centerline. The individual or sets of control nodes of a centerline may be freely translated, rotated or otherwise modified. The only requirement is that the second derivatives of the spline location parameters are recalculated after modification. This operation is very fast. The calculation of complicated geometries generally requires a high level of computational intensity or simplification. In the event-based models, the geometric construction (beyond the calculation of a centerline) is postponed to the end of the algorithm. This leads to fast calculation and manipulation of complicated geometric morphologies and associations represented as centerlines. Any properties may be attached to the centerline and interpolated along the entire length. In the fluvial event-based model, the channel width, local curvature, relative thalweg location, local azimuth and near bank velocity are included in the centerline. Other information including architectural element type, geometric parameters and additional property trends may be included. 3.1. Centerline complexes within event-based models A centerline complex is a grouping of genetically linked centerlines. Centerline complexes are characterized by their external and internal structure. The external structure is the interrelationship between centerline complexes. This may include stacking patterns, such as dispersive patterns in alluvial fans, vertical stacking in anastomosing reaches and nested channel belts in incised valleys. These patterns include important information with regard to the heterogeneity of a reservoir and should be included in fluvial models. The internal structure is the relation of centerlines within the centerline complex. This may include meander migration and resulting lateral accretion (LA) architectural elements or frequent avulsion and braiding. The internal structure is constrained by a set of centerline operations. 3.2. Centerline operations Centerline operations generate and modify centerlines in a manner that mimics fluvial processes. Architectural element models are calculated by sequentially applying the following centerline operations: (1) initialization, (2) avulsion, (3) aggradation, (4) migration and (5) neck cut off.

Fig. 1. Left—an illustration of a centerline association. Right—associated architectural element model. Models are continuous in distal direction, but are split to show cross section of architectural elements.

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

3.3. Centerline initialization

1673

4000. 1.8 3000. Y Coordinate (m)

The centerline initialization operator is applied to generate a centerline or to represent channel avulsion proximal of the model area. Given source and target locations and channel sinuosity, this operator initially generates a realistic centerline beginning at the source. A smooth translation is then applied to correct the location of the centerline to honor the target location (Fig. 2). A model is required to represent channel centerline morphology. The disturbed periodic model developed by Ferguson (1976) as a meandering river analog provides a realistic centerline model for a fluvial depositional environment.

1.7 1.6 1.5 1.4 1.3 1.2

2000.

1.1

1000.

2

fðsÞ þ

2h dfðsÞ 1 d fðsÞ ¼ ðsÞ þ 2 k ds ds2 k

(1)

where k is related to the primary wavelength k ¼ 2p/l, h the dampening factor (0oho1) and e(s) is the disturbance value. The physical analogue for this model is a pendulum dampened by air resistance and continuously hit by rocks. The discrete approximation of this model may be applied to efficiently calculate centerlines with intuitive parameters and with the periodicity and irregularity observed in fluvial channels. Sinuosity is related to the dampening factor, h, the wavelength through k, and the variance of the disturbance random variable, e. A decrease in dampening results in more regularity and higher sinuosity. Decrease in the disturbance variance increases the regular periodicity of the model. A series of channel centerlines are generated with varying sinuosity in Fig. 3. The interested reader is referred to Ferguson (1976) for the mathematical details of this model.

0. 0.

1000.

2000.

3000.

4000.

X Coordinate (m) Fig. 3. Example channel centerlines calculated with disturbed periodic model. From bottom to top sinuosity ¼ 1.1, 1.2, 1.3, y , 1.8.

3.5. Centerline aggradation Aggradation occurs when a channel deposits sediments in the channel and overbank environments. This process is represented by an incremental increase in the elevation of a centerline. The current implementation is to add a specified constant value to the elevation, z, parameter for all centerline control nodes. 3.6. Centerline migration

3.4. Centerline avulsion

S

T

Proximal

Transverse Coordinate (m)

Avulsion is the process of channel abandonment and the establishment of a new channel outside the active channel and is handled with the centerline avulsion operator. This operator creates a copy of a specific channel centerline and selects an avulsion node stochastically. A new downstream centerline segment is generated with same centerline sinuosity and the same geometric parameter distributions as the original centerline. The geometric parameters (e.g. channel width, levee height, etc.) of the new centerline are corrected, so that the properties are continuous at the avulsion location. The initial azimuth is specified as the azimuth of the tangent at the avulsion location. There is no constraint to prevent the avulsed centerline from crossing the original centerline distal of the avulsion location. The avulsion location is drawn from the distribution of centerline control nodes weighted by the curvature at each node. This rule integrates the concept that avulsion is more likely to occur at locations of high curvature.

Distal

The centerline migration operator is based on the bank retreat model. The application of the bank retreat model for realistic channel meander migration has been proposed by Howard (1992), applied to construct fluvial models by Sun et al. (1996) and extended to construct meandering fluvial models that approximately honor global proportions, vertical and horizontal trends by Lopez et al. (2001). In the ALLUVSIM algorithm, the meander migration along the centerline is standardized such that the maximum migration distance matches a user-specified value. This removes the significance of hydraulic parameters such as friction coefficient, scour factor and average flow rate, since only the relative near bank velocity along the centerline is significant. Hydraulic parameters are replaced by the maximum spacing of accretion surfaces.

4. Neck and chute cut offs Neck cut offs (an entire meander loop is abandoned) and chute channels (the channel cuts across the point bar) are generated with the neck and chute cut off operation. The control nodes are searched for nonadjacent control nodes that are within a channel width of each other. When this occurs, the segment between these control nodes is removed and coded as abandoned channel (Fig. 4 for an illustration of the cut offs).

5. Fluvial architectural elements Longitudinal Coordinate (m)

Fig. 2. An illustration of centerline initialization operator. A smooth translation is applied to an initial centerline (gray nodes) to correct target location, resulting in final spline fit (black nodes).

Once the centerlines have been created, they must be populated with architectural elements. The available architectural elements include overbank fines (FF), channel fill (CH), lateral

ARTICLE IN PRESS 1674

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

4000

Chute Cutoff Neck Cutoff

Y Coordinate (m)

3000

2000

1000

T ime Step 1

50

100

0 0

1000

2000

3000

4000

X Coordinate (m) Fig. 4. Migration of a channel centerline over 100 times steps (white to black). Source and target locations are kept constant. Note neck and chute cut offs and end point constraints.

LV_width

LV_width Ratio Cutbank:Pointbar

CH_width:depth

LVtop LV_height

CH_elevation

LV_depth

CH_depth

LVbase location_thalweg

CH base

Fig. 5. Channel and levee elements cross section and associated geometric parameters. Note LV geometry is eclipsed by channel geometry.

accretion, levee, crevasse splay and abandoned channel fill (FF(CH)), refer to the conceptual geometric model in Fig. 1. The geometries and associated parameters are discussed for each element. 5.1. Overbank fines (FF) elements FF elements are not modeled directly. The model space is initialized as FF elements. All subsequent architectural elements displace the background FF during model construction.

(24) and Fig. 10 in Deutsch and Wang, 1996). The CH geometry and associated parameterization are shown in Fig. 5. This geometry is consistent with the general form observed in meandering streams (Easterbrook, 1969). The CH element parameters; thalweg location, depth and width to depth ratio are calculated at control nodes along the centerline. Cubic splines are fit to these properties to allow for a smooth transition along the centerline and efficient interpolation at any centerline location. The CH element average width to depth ratio and average depth are drawn from Gaussian distributions specified in the input parameters. Conservation of cross-sectional area is honored along the centerline.

5.2. Channel fill elements The geometry of the CH element is parameterized by a centerline, relative thalweg location, stochastic depth and a width-to-depth ratio. The cross section channel geometry is based on the definition from Deutsch and Wang (see Eqs. (22)–

5.3. Lateral accretion elements The lateral accretion deposits are represented as the volume of channel abandoned during channel migration. The complicated

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

1675

geometries of LA elements have been demonstrated by Diaz-Molina (1993) and Willis (1993). LA elements are characterized by wedge channel fills distributed along the inside of meander bends.

L

W 5.4. Levee elements Significant reservoir quality facies may be represented by levees. The LV geometry and associated parameterization are shown in Fig. 5. The equations for the LVtop and base are shown in Eqs. (2) and (3) respectively.   d LV top ðd; sÞ ¼ LV height ðsÞ LV width ðsÞF   d þ CHelev ðsÞ (2)  exp LV width ðsÞF LV base ðd; sÞ ¼ CHelev ðsÞ  LV depth ðsÞ



LV width ðsÞF  d LV width ðsÞF



F cutbank ðsÞ ¼ 1:0 þ LV asym

cðsÞ jcmax j

cðsÞ jcmax j

l

(3)

where LVtop is the top elevation of the LV element and LVbase is the elevation of the base of the LV element at a location s along the centerline and location d orthogonal from the channel edge. CHelev(s), LVdepth(s), LVwidth(s) and LVheight(s) are the elevation of the channel and the depth, width and height of the LV element at the nearest location along the centerline, s. At locations with negative thickness (i.e., LVbase(s, d)4LVtop(s, d)) no LV element is placed. F is a scaling factor described below. The distribution of LV elements may not be uniform along the channel axis. Typically LV elements are more pronounced on the cut bank. This information is integrated with a LV width scaling factor, F, that scales the LV width to account for LV asymmetry on channel bends. The current implementation is applied to calculate the F multiplier with the following equations. F point bar ðsÞ ¼ 1:0  LV asym

w

CH

LV

CS

Streamline

Avulsion Node

Fig. 6. Crevasse splay architectural element geometric parameters. L is length of lobe, l is length along the centerline at which lobe has its maximum width, W is width along proximal edge.

Example, CS elements constructed with ALLUVSIM are shown in Fig. 7.

(4) 5.6. Abandoned channel (FF(CH)) elements (5)

where LVasym is a value between 0 and 1 that parameterizes the strength of the levee asymmetry, c(s) the local curvature along the centerline and cmax the maximum curvature along the centerline. LVasym value of 0 results in symmetric LV elements and a LVasym value of 1 results in no LV element on the point bar side and double width of LV element on the cut bank side at the location along the centerline with greatest curvature. 5.5. Crevasse splay (CS) elements CS elements may also have significant reservoir quality. For each centerline, the number of CS elements is drawn from a Gaussian distribution with mean and standard deviation supplied by the user. The location of each CS element along the centerline is drawn from a distribution of centerline locations, weighted by the curvature. Crevasse splays more likely occur at locations with high curvature, since high near bank velocities erode the confining LV elements. The CS elements are modeled as a series of lobes fit to low-sinuosity centerlines initiated from the crevasse location with initial azimuth normal to the channel centerline towards the cut bank. The number of lobes and the lobe parameters (Fig. 6) are drawn from Gaussian distributions with user supplied mean and standard deviation. Fewer large lobes reproduce sheet geometries, while many thin lobes reproduce inter fingering. This geometry is taken from the Lobesim algorithm by Deutsch and Tran (2002).

When channel abandonment is very abrupt (i.e., rapid avulsion and neck cut off) there is a strong contrast between the abandoned channel fills FF(CH) and CH elements, while slow abandonment leads to fining upward fills (Galloway and Hobday, 1997). In the current implementation FF(CH) elements form: (1) along cut off centerline segments, (2) in the channel reaches distal of avulsion locations and (3) in the last channel placed for a level. The user specifies the mean and standard deviation of the proportion of abandoned channels fill with fine-grained FF(CH) elements. For a proportion of 0 the abandoned channel is coded as CH element and for a proportion of 1 the entire abandoned channel is coded as FF(CH) element. For a proportion between 0 and 1, the abandoned channel fill is partitioned with proportional coordinates (refer back to the geometric model in Fig. 1).

6. Event schedule and conditional to areal and vertical trends The ALLUVSIM algorithm is able to reproduce a wide variety of reservoir styles. This algorithm may reproduce braided, avulsing, meandering channels and may reproduce geometries and interrelationships of a variety of reservoir types (PV, CB or sheet (SH) type) and fluvial style within any systems tract. The algorithm can be supplied with areal and vertical trends, distributions of geometric parameters, probabilities of events and architectural elements.

ARTICLE IN PRESS 1676

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

1000

Y Coordinate (m)

Y Coordinate (m)

1000

0

0 0

1000

CH 0

1000 X Coordinate (m)

X Coordinate (m)

LA LV

1000

1000 Y Coordinate (m)

Y Coordinate (m)

CS FF(CH) FF

0

0 0

1000 X Coordinate (m)

0

1000 X Coordinate (m)

Fig. 7. Example crevasse splay geometries calculated with ALLUVSIM algorithm: A—isolated inter-fingered elements, B—amalgamated lobe and inter-fingered elements, C—isolated lobe elements and D—amalgamated lobe elements.

6.1. Areal channel density trends The technique for honoring areal trends is to (1) construct a suite of candidate centerlines with the desired morphology, (2) superimpose each centerline on the areal trend model and calculate average relative quality and (3) for each centerline initialization draw from this distribution of candidate centerlines (without replacement) weighted by the average quality index. This technique is efficient since the construction of hundreds or thousands of centerlines is computationally fast. This technique is demonstrated in Fig. 8, 30 centerlines are drawn from a suite of 500 candidate centerlines for three different areal trends. Note the impact of the areal trend on the drawn centerlines.

6.2. Vertical channel density trends and aggradation schedule Vertical trends may be honored by constraining the aggradation schedule. The current implementation is to apply the trend within a user-defined number of constant elevation levels. Centerlines and associated architectural elements are generated at the lowest level until the NTG indicated by the vertical trend is reached for the model subset from the base of the model, z0, to the elevation of the first level, z1. Then the aggradation operator is applied to aggrade to the next level and the process is repeated through all user-defined levels. For the highest level, zn, the model is complete when the global NTG ratio is reached. The assignment of a few constant elevation levels to model aggradation is a simplification, although this technique does allow

for the flexibility to reproduce a variety of aggradation patterns. For example, few levels result in a stepped system with amalgamated reservoir elements separated by FF elements. The assignment of many levels and an overall low NTG results in isolated PV-type shoestring or CB-type labyrinth reservoirs. The assignment of many levels with an overall high NTG results in an amalgamated CB-type jigsaw reservoirs. With minor modification more complicated schedules may be applied that allow for continuous aggradation as opposed to discrete levels. Three example vertical trends with cross sections from the resulting models are shown in Fig. 9.

7. Conditioning event-based models to honor well data There are a variety of available methods that may be applied to condition complicated geologic models; (1) dynamically constrain model parameters during model construction to improve data match (Lopez et al., 2001), (2) iterative methods such simulated annealing and maximum a-posteriori selection (Deutsch, 1992), (3) pseudo-reverse modeling (Tetzlaff, 1990), (4) apply as a training image for multiple-point geostatistics (Strebelle, 2002), (5) direct fitting of geometries to data (Oliver, 2002; Viseur et al., 1998 and Shmaryan and Deutsch, 1998) and (6) convert the architecture into a continuous random function and apply traditional geostatistical methods for conditioning (Alapetite et al., 2005). Each of these techniques has significant limitations either in efficiency, robustness or the ability to retain complicated geometries and interrelationships. The event-based paradigm is amenable to a new method for the construction of conditional models. The procedure is: (1) code

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

1000

1677

1000

15 10 5

Y Coordinate (m)

Y Coordinate (m)

20

0 0

0 0

1000

0

X Coordinate (m)

1000 X Coordinate (m)

1000

1000

30 20 10

Y Coordinate (m)

Y Coordinate (m)

40

0 0

0 0

1000

0 X Coordinate (m)

X Coordinate (m)

1000

1000

1000

75 50 25

Y Coordinate (m)

Y Coordinate (m)

100

0 0

0 0

1000 X Coordinate (m)

0

1000 X Coordinate (m)

Fig. 8. Example areal trends in channel density and resulting centerlines. A and B—no areal trend supplied, C and D—a linear trend increasing in positive direction and E and F—a second-order trend increasing in positive direction. Note areal trend is a relative measure without units.

the well data as net or non-net, (2) construct the prior event-based model conditioned by all available soft information and constrained by rejection criteria to prevent significant contradictions with conditioning during construction, (3) sequentially update centerline complexes to honor identified net intervals and (4) apply image cleaning to remove small occurrences of well data and model mismatch and to honor specific architectural elements.

8. Net and non-net coding All architectural elements are grouped together as net and the background is coded as non-net (Fig. 10). This greatly simplifies the procedure for iteratively updating centerlines to honor well

data. The post-processing by image cleaning then corrects the model to honor specific architectural elements.

9. Gross well conditioning during model construction After the calculation of a new centerline, the associated architectural elements are calculated at all nearby well locations. The mismatch between the model and well data net elements and non-net elements is quantified by the number of cells correctly matched subtracted by the number of cells incorrectly matched. A threshold of acceptable mismatch is set by the user. If this threshold is exceeded then the centerline is rejected and a new centerline is calculated. A high threshold may speed model construction, but result in artifacts at well locations.

ARTICLE IN PRESS 1678

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

20 Z Coordinate (m)

Z Coordinate (m)

20.0

15.0

10.0

5.0 0

0.0 1.0

5.0

9.0

13.0

17.0

0

20.0

20 CH

Z Coordinate (m)

Z Coordinate (m)

1000 X Coordinate (m)

Relative Reservoir Quality

15.0

10.0

LA LV CS

5.0

FF(CH) FF

0

0.0 1.0

5.0

9.0

13.0

0

17.0

20.0

20 Z Coordinate (m)

Z Coordinate (m)

1000 X Coordinate (m)

Relative Reservoir Quality

15.0

10.0

5.0 0

0.0 1.0

5.0

9.0

13.0

17.0

Relative Reservoir Quality

0

1000 X Coordinate (m)

Fig. 9. Example vertical trends in channel density and resulting architectural element models. A and B—no areal trend supplied, C and D—a linear trend increasing in positive direction and E and F—a second-order trend increasing in positive direction. Source locations are drawn from a uniform distribution along proximal edge. Note vertical trend is a relative measure without units.

10. Updating centerline complexes to honor well data The model is updated by modifying centerline complexes to honor net element intercepts. For each net element interval identified in the well data, the nearest centerline association is found and updated in the following manner: (1) The horizontal position is corrected such that the net intercept thickness is within tolerance of the net interval thickness, (2) the vertical location is corrected such that the net intercept top matches the top of the net interval. Entire centerline complexes are corrected to preserve the relationships between centerlines within a centerline association. For example, if a centerline association includes a set of centerlines related by meander migration, the entire set of centerlines representing a point bar is shifted. If

individual centerlines were modified independently, this would change the nature of the centerline association.

10.1. Iterative procedure for updating centerline complexes A modification of centerline complexes has an impact on net element geometry. It would be difficult to directly calculate the precise translation of a centerline to result in the correct interval thickness at a well location. An iterative method is applied to correct the well intercept thickness to obtain approximate conditioning to well data. The thickness of the net from a centerline association is calculated at the vertical well location. If the thickness is less than indicated by the conditioning,

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

1679

Fig. 10. A cross section illustrating net element coding.

1000 Y Coordinate (m)

c 1 p

n1 n2

0 X Coordinate (m)

Distal

Δ(m)

Proximal

1 0

1000

p

c

n 1 ,n 2

s(m)

Fig. 11. An illustration of methods for updating centerline complexes to well data. For this example, there are two centerlines in the centerline association representing an avulsion event corrected to honor conditioning data (c). A—centerline is translated with a smooth difference distribution with the centerline held constant at control node of previous conditioning (p) and beyond (p41) and at last control nodes (n1 and n2). B—transverse correction with respect to location along the centerline.

then the centerline association is shifted towards the well location (Fig. 11). 10.2. Image cleaning for precise reproduction of architectural element intercepts The object-based framework of models generated with eventbased ALLUVSIM contain realistic fluvial features and have largescale curvilinear behaviors that are not easily created by other techniques; however, it is difficult to reproduce dense well control

with the placement of large-scale objects, and there are often slight mismatches at well locations. Consider the mismatches at the well locations in Fig. 12. It may be unrealistic and unnecessary to expect precise well conditioning, under the constraint of a very specific morphological model as shown by the schematic black line. Convergence to a precise solution may be prolonged or intractable, sacrificing algorithm robustness. Algorithm robustness may be improved by accepting a acceptable level of mismatch and then applying a smooth correction. The goal is to have geologic realism, as provided by ALLUVSIM, with perfect well data reproduction. Post-processing of the ALLUVSIM models to exactly reproduce well data will be accomplished with an image-cleaning algorithm, maximum a-posteriori selection postprocessing (MAPS-PP). For the original implementation of MAPS, the interested reader is referred to Deutsch (1998). Essentially, the original implementation randomly visits every location in a model, determines the most likely facies value based on the surrounding facies and reassigns the target location to the most likely facies (Fig. 12: in the case of dilation this would be the net facies and in the case of erosion this would be the non-net facies). The object-based models should only be changed in areas where the observed well data do not match the ALLUVSIM models. Cells that are candidates for a change in facies are identified as those within an ellipsoidal range from the cell–well mismatch. The cells at the mismatched locations will be visited first and a spiral search will be used until the ellipsoidal range is reached. The probability (p) of each facies prevailing at a particular cell location (u) is calculated based on a weighting function pðu; kÞ ¼

1 X wðu0 Þ  iðu; kÞ; S u0 2WðuÞ

k ¼ 1; . . . ; K

where S is a standardization constant, W(u) is a rectangular template of weights centered at the location under consideration and i(u; k) is the indicator of facies k at location u. There was much discussion on the weighting template in the original MAPS-PP paper. Based on the authors’ experience, a reasonably small template (5  5  5) works well.

ARTICLE IN PRESS 1680

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

Erosion

Dilation

0.8 + f

0.6 + f

Cell-well mismatch

0.4 0.2

Net Non-Net

Higher probability of non-net Replace with non-net facies

Cell-well mismatch

Higher probability of net Replace with net facies

Net Non-Net

Fig. 12. Handling cell-to-well data mismatch. A cross section of initial ALLUVSIM model is shown as gray (net) and white (non-net) cells and post processed using MAPS-PP. Well data is coded as net (think black line) or non-net (thin black line). Histograms indicate facies proportions calculated using MAPS-PP before adding distance correction f. Left—dilation is shown as 5 non-net cells will be recoded to net because of cell-to-well mismatch. Right—erosion is shown as 3 net cells will be recoded to non-net, because of cell-to-well mismatch.

Fig. 13. Before (above) and after (below) applying MAPS-PP for well data reproduction. Gray cells are considered net and black well data are considered non-net. Well data is exactly reproduced after applying MAPS-PP. Note that scale of figures would be dependent on site specific channel geometry and is not necessary in this schematic representation.

It is reasonable to conclude that the closer a cell is to a mismatch, the more likely that cell is to be changed to the well value. The probability that a cell is changed because of a mismatch should be one at the well location and decrease as the distance from the well increases. This is accomplished by adding a distance factor to the probability of changing a cell f ¼ ð1  dÞo where f is the amount to increase the probability of changing a cell, d is the standardized distance between 0 and 1, where d ¼ 0 at the well location and d ¼ 1 at the maximum distance from the cell. o controls how quickly this factor decreases with distance; a value of o ¼ 2 was found to be reasonable. f is only added to the probability if the cell under consideration is conflicting with the well data (i.e., has the opposite facies designation as the well data at the mismatch as in Fig. 12). The probability of the well facies prevailing at a location u within the search ellipsoidal range is p0 ðu; kÞ ¼ f þ pðu; kÞ The cell value is changed to the facies with the highest probability (Fig. 12). The algorithm is applied sequentially in a spiraling fashion away from the well. The original MAPS-PP algorithm was not

sequential; however, the goal here is for the erosion/dilation to be smooth away from the wells and not to allow unrealistic shortscale variations that would occur if a random path were used as per the original implementation. The algorithm is fast and wells are honored smoothly with few visual artifacts. Fig. 13 shows a typical result of cleaning an objectbased realization with the modified MAPS-PP algorithm.

11. Limitations in conditioning This technique entails the translation of large-scale elements to honor small scale data; therefore, it is only suitable for settings with sparse conditioning data. Settings with dense data will be intractable and may result in artifacts such as changes in architecture near the wells. The current implementation considers only net and non-net during model updating with post-processing applied to honor specific architectural elements. This may lead to some artifacts, although this feature may mimic remnant architectural elements. It is assumed that the wells are vertical. The implementation of deviated and horizontal wells would require modifications in the methods for loading, storing and upscaling well data to the model grid.

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

12. Program parameters The ALLUVSIM program is research quality code and is suitable for experimenting with the ideas presented in this paper. The code is from the Ph.D. thesis of Pyrcz (2004). The program follows GSLIB conventions for parameter interface, input and output files. An example parameter set is shown in Table 1 and a list of parameters with associated descriptions is provided below.

 

 Line 1: input file containing the channel fill element intervals      

from well data. The standard GSLIB/GEOEAS format is expected and the wells are assumed to be vertical. Line 2: columns in the well data file for the X, Y and Z (top and base) coordinates, the well number (used to identify different well intersections). Line 3: input file with relative horizontal trend in channel density. The file should be in GEOEAS format and GSLIB grid convention. Line 4: the column number for the horizontal channel density trend. Line5: input file with relative vertical trend in channel density. The file should be in GEOEAS format and GSLIB grid convention. Line 6: the column number for the vertical channel density trend. Line 7: the maximum number of centerlines, the constant elevation levels and a list of the associated levels. The former is applied to dimension centerline tables. The algorithm





1681

terminates when this number of centerlines are generated. The latter is applied to define the vertical spacing of channels. Line 8: the target net-to-gross ratio. The algorithm terminates when this net-to-gross ratio is exceeded. Line 9: the number of candidate centerlines calculated prior to model construction and the number of discretizations for spline interpolation between control nodes for locating the closest location along a centerline to a position in the model. Set the former several times larger than the maximum number of centerlines and set the latter at least 3, higher values will result in more precise architecture at the expense of computational time. Line 10: probability of avulsion proximal of the model (new centerline initialization) and the probability of avulsion within the model. One minus the sum of these probabilities is the probability of meander migration. Lines 11–15: channel fill architectural element geometric parameters. Parameterized by mean and standard deviation of a Gaussian distribution. J Line 11: primary azimuth of channel centerlines. Note current model assumes model proximal edge is x ¼ 0 and distal edge is x ¼ xmax. J Line 12: source location in Y coordinates. Source is located along the proximal edge of the model (e.g. X ¼ 0 for primary azimuth of 901). J Line 13: channel depth. J Line 14: width-to-depth ratio.

Table 1 An example parameter set and parameter descriptions for ALLUVSIM program. welldata.dat 1, 2, 3, 4, 7, 9 horitrend.dat 1 Verttrend.dat 1 90, 3, 2.0, 5.0, 7.0 0.3 100, 10 0.1, 0.1 90.0, 5.0 500.0, 50.0 5.0, 0.2 20.0, 2.0 1.3, 0.2 2.0, 0.1 160.0, 5.0 1.5, 0.2 0.3, 0.1 0.3, 0.1 2, 2 10, 2 300.0, 100.0 150.0, 30.0 50.0, 10.0 150.0, 30.0 0.000, 0.000 0.001, 0.005 0.5, 0.1 50.0, 20.0 200, 2.5, 5.0 200, 2.5, 5.0 50.0, 0.2, 0.4 13013, 0.05 alluvsim.out centerline.out

File with the CH element intercepts Columns for the well #, X, Y, Ztop and Zbase File with the horizontal trend Column with the horizontal trend File with the vertical trend Column with the horizontal trend Max # of events, # of levels and list of elevations NTG target # Candidate centerlines in look up table for horizontal trend matching, number of fine search intervals between centerline control nodes for placing architecture Probability of major avulsion proximal of model (independent centerlines) and probability of avulsion within the model. In the absence of avulsion, centerline meander is applied. CH element: azimuth* Source location*, mean distance along the proximal side and standard deviation ¼ 0 results in a point source Depth* from channel top to base at thalweg Width-to-ratio* Sinuosity* LV element: depth* from channel top to levee base Width* from channel centerline to levee margin Height* from channel top to levee crest Asymmetry factor*, symmetric levees—0 and double-sized cut bank levees—1 Thinning factor*, uniform cross section down dip—0 and thin to zero along centerline—1 CS Element: number of CS/centerline* Number of lobes/CS* Lobe length* Lobe maximum width* Lobe length to maximum length* Lobe proximal width* Lobe height-to-width ratio* Lobe depth-to-width ratio* FF(CH) element: proportion of channel fill* Maximum meander migration step* Grid: nx, xmin, xsize ny, ymin, ysize nz, zmin, zsize Random number seed and color increment on element codes File with output element model File with output centerlines (index, x, y and z)

Line Line Line Line Line Line Line Line Line

1 2 3 4 5 6 7 8 9

Line 10 Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line Line

11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36

Parameter file formatting conventions from second edition of GSLIB are followed, i.e., line number after ‘‘START’’ is significant and parameters are space delimited. These parameters result in a meandering model similar to Fig. 15.

ARTICLE IN PRESS 1682

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

Line 15: sinuosity. Based on a calibration of the Ferguson (1976) model. Lines 16–20: levee architectural element geometric parameters. Parameterized by mean and standard deviation of a Gaussian distribution. J Line 16: levee depth below top of channel fill. J Line 17: levee width from edge of channel fill. J Line 18: levee height above top of channel fill. J Line 19: factor for levee asymmetry on point bar and cut bank. For a value of 0, levees are symmetric and for a value of 1, levees are twice as wide on the cut bank side at the location of maximum curvature. J Line 20: factor for proximal to distal thinning along the centerline. For a value of 0, there is no thinning and for a value of 1 levee widths are doubled at the proximal edge and halved at the distal edge of the model. Lines 21–28: crevasse splay architectural element geometric parameters. Parameterized by mean and standard deviation of a Gaussian distribution. J Line 21: number of crevasse splays along a single channel centerline. J Line 22: number of lobes within a single crevasse splay. J





1000

Line 23: lobe length along the centerline. Line 24: lobe maximum width. J Line 25: lobe length along centerline to position of maximum length. J Line 26: lobe width at proximal edge. J Line 27: lobe height-to-width ratio. J Line 28: lobe depth-to-width ratio. Line 29: the geometric parameters for fine-grained abandoned channel fill element. The fraction of abandoned channel fill assigned as fine grained. Line 30: the maximum meander migration step. Translations calculated by the bank retreat model are standardized by this value. Line 31: the size of the model in the X-direction. Line 32: the size of the model in the Y-direction. Line 33: the size of the model in the Z-direction. Line 34: the random number seed (large odd integer) and the element code increment for differentiation of individual architectural elements. Line 35: output file containing the output-gridded architectural element realization. The realizations are written from the lower left corner and then realization-by-realization (X cycles fastest, then Y, Z and realization number). J J

      

Y Coordinate (m)

Y Coordinate (m)

1000

CH

LA 0

0 0

1000

0

X Coordinate (m)

1000

LV

X Coordinate (m) CS 20

FF(CH) FF

Y Coordinate (m)

Z (m)

1000

0 0

X Coordinate (m)

1000

0 0

1000 X Coordinate (m)

Fig. 14. An example ALLUVSIM CB-type jigsaw reservoir model. A—plan section (z ¼ 5 m), B—plan section (z ¼ 10 m), C—all centerlines (grey scale from 1 ¼ white to n ¼ black) and D—cross section (x ¼ 10 m). Note grey scale assignment for architectural elements is varied to aid in differentiating amalgamated elements. Note meandering and braided features in centerlines.

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

 Line 36: output file containing the centerlines applied to

other conditioning methods and (3) model heterogeneities within architectural elements. The ALLUVSIM algorithm is tailored to the fluvial depositional settings. The centerline planforms, centerline operations and attached architectural elements are based on this setting. This modeling methodology is well suited to depositional settings that inherent their heterogeneities from flow events. A natural extension would be to apply this methodology to deepwater settings (Pyrcz et al., 2005; Pyrcz and Strebelle, 2006, 2008). The centerline meander operation may be modified to constrain the degree of swing and sweep and meander step size may be set large. This would allow for the modeling of stepped meander migration with limited sweep as observed in deepwater channels. The lobe architectural elements applied in CS elements may be applied to construct frontal splays, and unconfined sheets and outer levees may be added. With pixel-based simulation approaches, there is no genetic linkage between model cells. Object-based methods link model cells into objects and may provide some linkage between objects. The event-based approach constructs models that have fully linked architecture; elements, complexes and complex sets. A channel fill may transition into lateral accretion, levee, crevasse splays and overbank from center to periphery of a centerline association. Continuous properties may be informed by hierarchical trend models based on this advanced linkage (Pyrcz et al., 2005) or further refinement of facies categories, such as channel lags and mud drapes may be included (Pyrcz and Deutsch, 2004). Allogeneic controls may be introduced with trends in parameters to produce complexes and complex sets.

construct the architectural element model.

13. Example event-based fluvial models The jigsaw reservoir type is characterized by high net-to-gross with no major gaps (Galloway and Hobday, 1997 and Miall, 1996). These reservoirs form from coarse-grained meandering and braided fluvial styles. In these coarse-grained systems, the channels typically have high width-to-depth ratios. An example jigsaw reservoir model is shown in Fig. 14. Note the braided and meandering centerline complexes and the formation of FF(CH) baffles. In Fig. 15, a SH-type model based on a distal setting is shown. The centerline complexes show a high degree of meandering and little avulsion. A prior and updated model is shown in Fig. 16. The centerlines include braided low-to-high sinuosity morphology. A single well is included. The morphology of the centerlines is preserved while the well intercepts are honored.

14. Future work The ALLUVSIM algorithm provides the basic building blocks to fluvial event-based simulation. These tools may be rapidly extended to (1) model different depositional settings, (2) integrate

1000

1683

Y Coordinate (m)

Y Coordinate (m)

1000

CH LA

0

0 0

1000

0

X Coordinate (m)

1000 X Coordinate (m)

CS

20

FF(CH) FF

Y Coordinate (m)

Z (m)

1000

LV

0 0

1000 X Coordinate (m)

0 0

1000 X Coordinate (m)

Fig. 15. An example ALLUVSIM distal SH-type reservoir model. A—plan section (z ¼ 5 m), B—plan section (z ¼ 10 m), C—all centerlines (grey scale from 1 ¼ white to n ¼ black) and D—cross section (x ¼ 10 m). Note grey scale assignment for architectural elements is varied to aid in differentiating amalgamated elements.

ARTICLE IN PRESS 1684

M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

Well 2

18.0

Z Coordinate (m)

Z Coordinate (m)

Well

13.0 8.0

18.0 CH

13.0 8.0

LA 3.0

3.0

-2.0

-2.0 0.

0.

200. 400. 600. 800. 1000.

1000.

1000.

800.

800.

600. Well 400.

200. 400. 600. 800. 1000. Y Coordinate (m)

Y Coordinate (m)

Y Coordinate (m)

Y Coordinate (m)

CH 600. Well LA

400. 200.

200.

0.

0. 0.

200.

400.

600.

800. 1000.

X Coordinate (m)

0.

200.

400.

600.

800. 1000.

X Coordinate (m)

Fig. 16. An example conditional event-based model using ALLUVSIM with CH and LA elements. (a) cross section before conditioning, (b) cross section after conditioning, (c) plan view before conditioning and (d) plan view after conditioning showing movement of centerlines around the well.

15. Conclusions The event-based approach is a flexible and efficient tool for the construction of stochastic fluvial models. The building block approach allows for the modeling of a variety of fluvial reservoir styles, including the complicated architectures of CB-type fluvial reservoirs. Event-based models may be calculated based on all available soft geologic information and then updated to honor hard well data. ALLUVSIM provides a useful public domain research code for geologic modelers considering event-based architectural element modeling. References Alapetite, A., Leflon, B., Gringarten E., Mallet, J., 2005. Stochastic modeling of fluvial reservoirs: the YACS approach. In: Proceedings Society of Petroleum Engineers Annual Technical Conference and Exhibition, Society of Petroleum Engineers, Dallas, USA. SPE Paper 97271. Cojan, I., Fouche, O., Lopez, S., 2004. Process-based reservoir modelling in the example meandering channel. In: Leuangthong, O., Deutsch, C. (Eds.), Geostatistics Banff 2004. Springer, Netherlands, pp. 611–619. Collinson, J.D., 1996. Alluvial sediments. In: Reading, H.G. (Ed.), Sedimentary Environments: Processes, Facies and Stratigraphy. Blackwell Science, London, pp. 37–82. Deutsch, C.V., 1992. Annealing techniques applied to reservoir modeling and the integration of geological and engineering (well test) data. Unpublished Ph.D. Dissertation, Stanford University, Stanford, CA, p. 306. Deutsch, C.V., Wang, L., 1996. Hierarchical object-based stochastic modeling of fluvial reservoirs. Mathematical Geology 28 (7), 857–880. Deutsch, C.V., 1998. Cleaning categorical variable (lithofacies) realizations with maximum a-posteriori selection. Computers & Geosciences 24, 551–562.

Deutsch, C.V., Tran, T., 2002. FLUVSIM: a program for object-based stochastic modeling of fluvial depositional systems. Computers & Geosciences 28, 525–535. Diaz-Molina, M., 1993. Geometry of lateral accretion in meander loops: examples from the Upper Oligocene-Lower Miocene, Loranca Basin, Spain. In: Marzo, M., Puigdefabregas, C. (Eds.), Alluvial Sedimentation-Special Publication #17, International Association of Sedimentology. Blackwell Science, London, pp. 115–131. Easterbrook, D.J., 1969. Principles of Geomorphology. McGraw-Hill, New York, p. 462. Ferguson, R.I., 1976. Disturbed periodic model for river meanders. Earth Surface Processes 1 (4), 337–347. Galloway, W.E., Hobday, D.K., 1996. Terrigenous Clastic Depositional Systems. Springer, New York, p. 489. Howard, A.D., 1992. Modeling channel migration and floodplain sedimentation in meandering streams. In: Carling, P.A., Petts, G.E. (Eds.), Lowland Floodplain Rivers. Wiley, New York, pp. 1–42. Jin, Y., Liu, D., Luo, C., 1985. Development of Daqing oil field by waterflooding. Journal of Petroleum Technology 37 (2), 269–274. Lopez, S., Galli, A., Cojan, I., 2001. Fluvial meandering channelized reservoirs: a stochastic and process-based approach. In: Proceedings Annual Conference of the International Association of Mathematical Geologists. International Association of Mathematical Geologists, Cancun, Mexico, CD-ROM. Miall, A.D., 1996. The Geology of Fluvial Deposits. Springer, New York, p. 582. Mossop, G.D., Flach, P.D., 1983. Deep channel sedimentation in the lower cretaceous McMurray formation. Sedimentology 30, 493–509. Oliver, D.S., 2002. Conditioning channel meanders to well observations. Mathematical Geology 34, 185–201. Patterson, P.E., Jones, T.A., Donofrio, C.J., Donovan, A.D., Ottmann, J.D., 2002. Geologic modelling of external and internal reservoir architecture of fluvial depositional systems. In: Amrstrong, M., Bettini, C., Champigny, N., Galli, A., Remarce, A. (Eds.), Geostatistical Rio 2000. Kluwer Academic Publishing, Norwell, pp. 41–52. Pyrcz, M.J., 2004. Integration of geologic information into geostatistical models. Unpublished Ph. D. Dissertation, University of Alberta, Edmonton, Canada, p. 296.

ARTICLE IN PRESS M.J. Pyrcz et al. / Computers & Geosciences 35 (2009) 1671–1685

Pyrcz, M.J., Deutsch, C.V., 2004. Stochastic modeling of inclined heterolithic stratification with the bank retreat model. In: Proceedings of the 2004 Canadian Society of Petroleum Geologists. Canadian Well Logging Society and Canadian Heavy Oil Association Joint Convention (ICE2004), Calgary, Canada. Pyrcz, M.J., Deutsch, C.V., 2005. Conditioning event-based fluvial models. In: Leuangthong, O., Deutsch, C. (Eds.), Geostatistics Banff 2004. Springer, Dordrecht, The Netherlands, pp. 135–144. Pyrcz, M.J., Catuneanu, O., Deutsch, C.V., 2005. Stochastic surface-based modeling of turbidite lobes. American Association of Petroleum Geologists Bulletin 89 (2), 177–192. Pyrcz, M.J., Strebelle, S., 2006. Event-based geostatistical modeling of deepwater systems. In: Slatt, R.M., Rosen, N.C., Bowman, M., Castagna, J., Good, T., Loucks, R., Latimer, R., Scheihing, M., Smith, R. (Eds.), In: Proceedings of the Reservoir Characterization: Integrating Technology and Business Practices: Gulf Coast Section—The Society for Sedimentary Research Twenty-Sixth Annual Research Conference, Houston, USA, pp. 893–922. Pyrcz, M.J., Strebelle, S., 2008. Event-based geostatistical modeling. In: Ortiz, J., Emery, X. (Eds.), Geostatistics Santiago 2008. Springer, The Netherlands, pp. 1143–1148. Shmaryan, L., Deutsch, C.V., 1999. Object-based modeling of fluvial/deepwater reservoirs with fast data conditioning: methodology and case studies. In: Proceedings of the Society of Petroleum Engineers Annual Technical Conference and Exhibition, Society of Petroleum Engineers, Houston, USA, SPE Paper 56821.

1685

Strebelle, S., 2002. Conditional simulation of complex geological structures using multiple-point statistics. Mathematical Geology 32 (9), 2937–2954. Sun, T., Meakin, P., Josang, T., 1996. A simulation model for meandering rivers. Water Resources Research 32 (9), 2937–2954. Tetzlaff, D.M., 1990. Limits to the predictive ability of dynamic models the simulation clastic sedimentation. In: Cross, T.A. (Ed.), Quantitative Dynamic Stratigraphy. Prentice-Hall, Oxford, pp. 55–66. Thomas, R.G., Smith, D.G., Wood, J.M., Visser, J., Calverley-Range, A., Koster, E.H., 1987. Inclined heterolithic stratification—terminology, description, interpretation and significance. Sedimentary Geology 53, 123–179. Viseur, S., Shtuka, A., Mallet J.L., 1998. New fast, stochastic, boolean simulation of fluvial deposits. In: Proceedings of the Society of Petroleum Engineers Annual Technical Conference and Exhibition, Society of Petroleum Engineers, New Orleans, USA, pp.110–121. Wietzerbin, L.J., Mallet J.L., 1993. Parameterization of complex 3D heterogeneities: a new CAD approach. In: Proceedings of the Society of Petroleum Engineers Annual Technical Conference and Exhibition, Society of Petroleum Engineers, Houston, USA, SPE Paper 26423. Willis, B.J., 1993. Interpretation of bedding geometry with ancient point-bar deposits. In: Marzo, M., Puigdefabregas, C. (Eds.), Alluvial Sedimentation— Special Publication Number 17, International Association of Sedimentologists. Blackwell Scientific, Oxford, pp. 101–114.