- Email: [email protected]
Building 2-D stratigraphic and structure models from well log data and control horizons
Computers & Geosciences 29 (2003) 557–567
Building 2-D stratigraphic and structure models from well log data and control horizons Maulin D. Patel, George A. McMechan* Center for Lithospheric Studies, University of Texas at Dallas, P.O. Box 830688, 2601 N. Floyd Road, Richardson, TX 75083-0688, USA Received 25 April 2002; received in revised form 3 December 2002; accepted 15 December 2002
Abstract An algorithm to build a gridded 2-D seismic velocity (or any other physical property) model from well log data and control horizons is developed. Interpolation of well log data onto a 2-D grid uses inverse distance weighting or linear interpolation, guided by the shape of the control horizons that are predefined from seismic or other 2-D constraints. A key feature of the models is that they may contain layers that are truncated by unconformities or at faults, or that lap out smoothly at their tops and bottoms. Abrupt or smooth terminations are controlled by user flags. Applications are illustrated using resistivity and acoustic impedance log data from the Blake Ridge (offshore Carolina) and using a complex structure produced by overthrust tectonics in western Canada. Geologically reasonable models can be produced only if there are sufficient wells to sample every salient element in the model and sufficient control horizons to define the lateral character of the structures at the required level of detail. r 2003 Elsevier Science Ltd. All rights reserved. Keywords: Model building; Stratigraphic interpolation; Well log interpolation; Structural control
1. Introduction Time-to-depth conversion, numerical modeling, and depth imaging of seismic data require geologically and geophysically realistic velocity models (Berkhout and Rietveld, 1994; Rietveld and Berkhout, 1994). Increasing quality of seismic data and survey accuracy demand correspondingly high-quality velocity models. Models can be parameterized and interpolated in a variety of ways. Tomographic models are usually either constant property layers (Cai and McMechan, 1999) or regular grids (Bishop et al., 1985); recent advances include irregular grids that adapt with finer sampling as data resolution increases (Sambridge et al., 1995; Vesnaver, 1996; Curtis and Snieder, 1997; Weber, 2001). Withinlayer gradients are also possible (Kabir and Verschuur, 2000). A number of sophisticated commercial 3-D model builders also exist (Schultz and Canales, 1997; Mallet, *Corresponding author. Tel.: +972-883-2424; fax: +1-972883-2829. E-mail address: [email protected] (G.A. McMechan).
1992a, 2002; F.alt et al., 1991). Most of these are directed toward construction of either macro models for input to imaging (Clar et al., 1996), or detailed petrophysical models as part of seismic inversion (Bee and Jacobson, 1984), or for input to reservoir simulations (Aziz and Settari, 1979; Gundeso. and Egeland, 1990). Other applications include modeling stratigraphic channels and detecting faults (Mallet, 2002). Model building explicitly requires interpolation to provide consistent parameter values between locations where there are measurements. In the past, this has usually been done with statistical models and various forms of kriging (e.g. Hohn, 1999; Xu et al., 1992; Deutsch and Journel, 1992; Mallet, 2002). However, where deterministic constraints are available (for example, from seismic or ground penetrating radar data), it is possible to produce models with deterministic geometry. Examples of comparisons of statistical and deterministic models are given by Torres-Verdin et al. (1999) and Helgelsen et al. (2000). Where only widely spaced well logs are available, it is impossible to adequately characterize lateral changes in lithology because of
0098-3004/03/$ - see front matter r 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0098-3004(03)00039-6
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
558
insufficient sampling. This limits the features that can be reliably interpolated to those of the largest scales. Where fairly continuous profiles are available (such as seismic data) the observed structural geometry can be inserted as constraints during interpolation of the well log data. Many deterministic interpolations (including for irregular grids) already exist (e.g., Akima, 1978; Sibson, 1981; Mallet, 1992b). These tend to be rather complicated, and the commercial versions are very expensive. Thus, it seems worth to extract only the essential concepts and implement them in a small, flexible, easyto-use format. This paper describes a procedure for building 2-D velocity (or other physical property) models, given well log data and control horizon geometries. The control horizons are typically obtained from migrated seismic depth sections. Measurements of some physical property (say seismic velocity, density, resistivity or porosity) constitute the well log data. Well data provide depths to key formation boundaries from drilling data or from well logs. The algorithm interpolates well log data so as to find reasonable values of these properties at all points in a 2-D grid by following the spatial geometry of the control horizons. Linear interpolation or inverse distance weighting are used for interpolation. Linear interpolation uses data only from the wells to either side of a specified location; inverse distance weighting involves contributions from all wells. A variety of layer truncation types are included, and the algorithms can be extended to 3-D.
2. Model description A vertical earth section is represented on a 2-D plane (Fig. 1). The control horizons are lines that go from one side of the model to the other to define the layer shapes; these will generally be a combination of sedimentologi-
X
cal bounding surfaces and other internal features. Vertical lines are the positions of wells that provide the well log data to be interpolated. Interpolation is between the well logs, using ties between the wells specified by the predefined control horizons and flags that specify how the interpolation is to be done. The final output is a 2-D grid of interpolated values. 2.1. Model specifications The number and positions of the vertical model control lines are arbitrary, but they should be placed to provide accurate descriptions of the model geometry. The X positions of the vertical control lines include the locations of all the wells. The algorithm has an option of omitting any well(s) from the model interpolation. The model can extend beyond the leftmost and rightmost X positions of the vertical control lines. The depth of each control horizon is specified on each vertical control line. The uppermost control horizon (representing the air– earth interface) can be flat or curved and is not necessarily at depth ¼ 0: The input well data points start from control horizon # 1 (not necessarily from depth ¼ 0). The depth extents of the well data are arbitrary. The output grid spacing in the vertical direction (DZ) is equal to the vertical spacing between samples in the well log data, and is the same in all well logs. All control horizons extend from left to right across the entire section. Truncated model features are defined by curves that also cross the entire section, but coincide with other control horizons beyond points where the model features are truncated. The control horizons are numbered from top to bottom in increasing order (Fig. 1). Control horizons are single-valued in depth and do not cross each other. Effectively multivalued geological horizons are composed by connecting
Position Well #1
Well #2
Well #3
Depth
Z
Control Horizon 4
Control Horizons
Control Horizon 5
Fig. 1. A vertical earth section on a 2-D plane with horizontal coordinate X and vertical coordinate Z: Control horizons and well logs provide geological constraints and geophysical data to be interpolated, respectively.
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
segments of single-valued horizons. The geometry of the control horizons conform to the geological environment and may, for example, be provided by seismic reflection profiles.
3. The algorithm The parameters input to the algorithm are: the number of output grid points in the X (horizontal) and Z (vertical) direction, the maximum number of samples in the well log data among all the wells at a depth increment equal to the vertical grid spacing (DZ), the grid spacing in the horizontal direction (DX ), the grid spacing in the vertical direction (DZ), the number and positions of the well(s) in the X direction, the amount of smoothing to apply to the well logs, the number and coordinates of vertical control lines, the boundary treatment parameters for all the control horizons (described below), and the well log data for all the wells. The task is to interpolate the well log data to generate consistent values of the measured properties at all points in a 2-D grid by following the geometry defined by the control horizons. Linear interpolation or inverse distance weighting can be chosen for the interpolation. The output is a 2-D gridded distribution of the estimated physical parameter on a vertical section through the earth. Step-by-step execution of the algorithm is explained below. 3.1. Outline of the algorithm The first step is to compute the depths of the control horizons at each of the vertical control lines, which include (but are not limited to) all of the wells. When the X coordinate of a control point coincides with a well position, the Z coordinates constitute the control horizon depths at that well. Where a geological feature has pinched out or has been truncated, the corresponding control horizons are superimposed on other control horizons that do extend completely across the model; thus, two or more control horizons may be at the same depth at a well, but are not allowed to cross. In Fig. 1, wells #1 and #3 do not have grid points between control horizons 4 and 5, so if well #2 is omitted from consideration, interpolation cannot be performed between horizons 4 and 5. If the number of grid points between a pair of control horizons at a particular well is zero, then the corresponding well data are assigned the value 1: When the control horizon #1 is not at depth ¼ 0; grid points are inserted above first control horizon. These points will not have any data values from wells and so are also assigned the value of 1: Similarly, when a well does not extend to the grid bottom, data points inserted below the bottom of the well are also assigned the value 1: Thus, the processed well data cover the
559
entire depth of the model, whereas the actual well data may not. When a data value of 1 is found during the interpolation, that data value is ignored as if that well was not present. The algorithm works from top to bottom and from left to right through all the grid points. For each horizontal position of the grid point, there are four possibilities: the grid point is to the left of the leftmost well position, or the grid point is to the right of the rightmost well position, or it superimposes on a well position, or it is between two well positions. For inverse distance weighting the first, second and fourth possibilities are treated in the same way and an inverse distance weighting routine is called. Where a grid point superimposes on a well position, the data values used are those of the well data at that location (no interpolation). When linear interpolation is specified, a linear interpolation routine is called. All the grid points to the left of the leftmost well are assigned the data values of that well (no interpolation). Similarly all the grid points to the right of the rightmost well are assigned the well data values at the same stratigraphic level in that well (no interpolation). For the third and fourth possibilities, the linear interpolation routine finds the interpolated data values. To interpolate a value at point a (Fig. 2A), the values at point b in well #1 and point c in well #2 are used; a–c are all at the same relative vertical position between control horizons j 1 and j: For example, if a is 0.3 of the increment from j 1 to j; then we can find the data values at b and c; which are also 0.3 of the increment from j 1 to j; at wells #1 and #2. This is done by subsampling the well data so there are the same number of points between the two given horizons at all wells. Layer truncations are a special case that is described in Section 3.4. 3.2. Linear interpolation Consider linear interpolation for the representative point a in Fig. 2A. The value assigned at point a is given by Va ¼
Da Db ðVc Vb Þ þ Vb ; Dc Db
where Va ; Vb and Vc are the data values at points a; b; and c; respectively. Da ; Db ; and Dc are the distances of points a; b and c from the origin, respectively. 3.3. Inverse distance weighted interpolation Consider the representative point d on the dashed interpolation trajectory in Fig. 2B. The value assigned to the point d is given by P Wi Vi Vd ¼ Pi ; i Wi
560
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 Position Control horizon j-1
Depth
Well #1
c a
b
Control horizon j
(A) Well #3
Depth
Well #2
Well #4
Well #5
C ontrol Region 1 e
g Region 2
f
d
(B) Fig. 2. Geometry for (A) linear interpolation and (B) weighted interpolation. In (A) dashed line maintains same relative position between control horizons j 1 and j: In (B) wells #3 and #5 have no data values for region 1, so for all grid points in region 1, only data values of well #4 are considered in interpolation. Interpolation at point d; on dashed interpolation trajectory efg, uses data values at f and g for linear interpolation, or at e–g for inverse weighted interpolation. Interpolations of grid points to left of well #3 and to right of well #5 use data values from only well #3 and well #5, respectively, for linear interpolation, or weighted values from all wells, for inverse weighted interpolation.
where the sum is over the contributions from each of the available wells (e.g. i ¼ 3; 4; 5). Vi is the data value in the well i at the same relative position that point d is between the control horizons above and below d; and Wi is the weight given to Vi : The sum in the denominator normalizes the weights so that their sum is 1.0: Wi ¼ ðDid Þp ; where Did is the distance from well i to point d; and p is the user-selected power weighting of the contributions. To avoid division by zero, the observed data values are used at the wells (where Did ¼ 0). A data value of 1 is ignored as if that well was not present in the 2-D grid, so the weight of that well is not calculated and it will have no impact on the weighted value. 3.4. Truncated control horizons All the control horizons extend from left to right across the entire section. To handle the problem of smaller curve segments, which do not cross the entire section, they are superimposed on control horizons that do extend from left to right. For example, in Fig. 3A, horizon db ends at b; but it needs to be specified as dbgc. Similarly, horizon eg ends at g; but is specified as efgc. With the interpolation described above, if a control horizon is truncated, then the interpolated trajectories will converge at a point. For example, in Fig. 3A, the dashed interpolated trajectories between the input
control horizons dbgc and efgc will all converge at point g: If bg is a stratigraphic pinchout that occurred during the initial deposition of the layer, then this convergent interpolation geometry is geologically reasonable. However, it is also necessary to consider the alternative that horizon abc is an erosional unconformity that has truncated horizons db and eg; in this case, the appropriate geometry is obtained by extrapolation past b; guided by the single existing horizon fg (Fig. 3B). The data values along the vertical line segment bf are replicated vertically at each horizontal position along the curve fg (see the dashed extrapolation trajectories in Fig. 3B; values extrapolated above bc are discarded). Layers that are truncated at their bottom end, for example by a fault, can be similarly extrapolated to produce the desired corresponding geometry. Thus, there are three possible behaviors that the user can specify: layer convergence, layer truncation at its upper end, or layer truncation at its bottom end. The decision on which choice is appropriate should come from a priori information from the geological environment. A similar, but more complicated solution to the layer truncation problem is illustrated by Mallet (2002); this involves one real and one virtual control horizon.
4. Examples In this section we present a series of examples from both synthetic and field data that illustrate the
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 Position
561
Position h
Depth
a
c
b
A
c
b g
g d
D f
f e
(A)
E
(B)
Fig. 3. Handling of truncated control horizons. (A) With layer convergence, dashed interpolated trajectories between input control horizons dbgc and efgc will converge at point g: (B) With truncation, appropriate geometry is obtained by extrapolation, guided by single existing horizon fg.
capabilities of the software. The examples show the behavior of many of the options described above. We show interpolation of resistivity and acoustic impedance well logs from the Blake Ridge, and compare with the results from a seismic inversion. This region is stratigraphically complicated, but structurally simple. Then, we apply the software to building a model from a structurally complex overthrust region in western Canada.
4.1. Synthetic example This example uses synthetic data and control horizons that illustrate the range of capabilities of the software. Thus, this model is not constrained to being geologically reasonable. Fig. 4A shows the locations of 16 vertical control lines (three of which are wells), 24 control horizons, and the range of data values at the wells. The boundary treatment parameters for control horizons and the interpolation algorithms to be used are specified by the user. Fig. 4B shows the estimates obtained by linear interpolation. Note the ability to deal with layer convergence (at a), truncation of a layer at its upper end (at b), and layer truncation at its bottom end (at c). Weighted interpolation (Fig. 4C) provides estimates when the X positions of the control points do not extend across the entire model. Removing the first and last vertical control lines results in the control horizons being extended horizontally to the model edges using the same Z positions as at the leftmost and rightmost input vertical control lines, respectively. Fig. 4D shows the effect of changing boundary treatment parameters from truncation to smooth convergence everywhere, and omitting well #2 from consideration. Omission of well #2 results in region d (which gets data only from well #2) turning white (i.e. becomes undefined), and specifying layer convergence of all boundaries results in all interpolation trajectories converging at points. Thus, knowing either layer shapes or having logs of the values to be interpolated is not sufficient; both must be present.
The interpolation procedures described above cannot by themselves produce geologically valid models. The models are only as complete as the constraints that are provided. For example, if only two well logs that are 5 km apart are given, the interpolation necessarily produces continuity across B2:5 km (if the quantity being interpolated has different values at the corresponding relative depths in the wells) or across all 5 km (if the quantity log interpolated has similar values at the corresponding relative depths). Additional lateral detail can be inserted only by providing additional constraints, for example, on the geometry (perhaps from seismic control) between the holes (e.g. Fig. 4) or by requiring the interpolation to confirm to a statistical model that defines a distribution of expected lateral dimensions. Thus, this approach has the advantage that the result is faithful to the type and resolution of the information that is input.
4.2. Stratigraphic examples The examples in this section are based on well logs from Ocean Drilling Program (ODP) Leg 164 holes 994D, 995B, and 997B on the Blake Ridge off the east coast of North America. Details of the data acquisition and processing are presented by Paull et al. (1996), Dillon et al. (1996), and Lu and McMechan (2002). Fig. 5 shows the locations of the three holes, the five control horizons that were extracted (as coherent reflections) from a seismic line through the holes, and the resistivity well log profiles measured in the holes. With these as input, we construct the interpolated sections in Fig. 6 using four different interpolations. The main differences are near the edges of the section. Linear interpolation (Fig. 6A) includes extrapolation, before the first and beyond the last holes, with the same data values as at those holes and so provides estimates only over the stratigraphic depths at which the logs exist (compare with Fig. 5). Weighted interpolations (Figs. 6B–D) provide estimates over all stratigraphic
Fig. 4. Synthetic example: (A) Locations of wells, control horizons, vertical control lines and data values at wells. Numbers x þ y þ z represent velocities at that well between a pair of control horizons in order of increasing depth. (B) Result of linear interpolation. (C) Result of inverse distance weighting with exponent value of 2 and omitting first and last vertical control lines. (D) Result of inverse distance weighting with exponent value of 2 and omitting well #2; specifying convergent layer boundaries results in layer convergence everywhere.
562 M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
563
Fig. 5. Stratigraphic example. Figure shows locations of wells, control horizons and resistivity log profiles measured in wells of ODP Leg 164 holes 994D, 995B, and 997B on Blake Ridge off east coast of North America. Corresponding interpolated sections are in Fig. 6.
depths at which a log exists in any hole. Away from the holes, the weighted extrapolations approach the average value of the logs, at the corresponding stratigraphic depth, at a rate that depends on the chosen exponent; for example, at the right side of Fig. 6D (with exponent 2), the rate of change is more rapid that in the same region of Fig. 6B (with exponent 0). Figs. 7A–C show the estimates obtained by linear and weighted interpolations of the acoustic impedance (the product of P-wave velocity and density) logs using the same control horizons as used in Fig. 6. For comparison, Fig. 7D shows the acoustic impedance estimated by inversion of the seismic section between the holes (after Lu and McMechan, 2002). The overall behaviors of the various estimates are similar; in particular, the low-impedance layer near 550–600 m depth (which also has high resistivity in Fig. 6) is clearly visible in both the estimates obtained by interpolation and by seismic inversion. This layer was interpreted by Lu and McMechan (2002) as containing free methane. 4.3. Structural example In addition to stratigraphic interpolation/extrapolation, as illustrated in the foregoing subsection, the software has utility in building structurally complicated (multivalued) models. Fig. 8A contains the vertical control positions and the layer shapes for building a velocity model for a data set used in SEG workshops on structural imaging in 1995 and 1996. Results based on this model have been presented by many authors, including Bevc (1997), Wu et al. (1998), and Zhu and Lines (1998). Fig. 8B contains the output 2-D velocity model. Note the ability to deal with multivalued
horizons (e.g., in region a of Fig. 8B), including faults (e.g., in region b of Fig. 8B), layer truncations (e.g., in region c of Fig. 8B), and pinchouts (e.g., in region d of Fig. 8B). Thus, all the main types of structural relationships can be included. Although Fig. 8A appears complicated, its generation is facilitated by a graphical user interface that allows points to be added, moved or deleted on a computer monitor, with mouse commands. This makes use in subsequent applications (e.g., updating velocities or layer geometries in iterative prestack migration) not only feasible, but also user-friendly. Extension to 3-D model building is also feasible, requiring mainly that the control horizons be specified as 3-D surfaces rather than as 2-D lines.
5. Summary We have developed and illustrated the operation of an algorithm for building 2-D models by interpolation and extrapolation of well log data. The interpolations and extrapolations are constrained by control horizons that are defined by seismic (or other) data that provide geometrical guides between the holes. Interpolation can be linear (between holes) or inverse distance weighting from all holes. Examples include interpolation between logs from a stratigraphically complex, but structurally simple, region off the east coast of North America, and building a structurally complex model for an overthrust zone in eastern Canada. Geologically reasonable models can be produced only if there are sufficient wells to sample every salient element in the model and sufficient
Fig. 6. Comparison of resistivity estimates: (A) result of linear interpolation, (B) result of inverse distance weighting with exponent value of 0, (C) result of inverse distance weighting with exponent value of 1, and (D) result of inverse distance weighting with exponent value of 2. White portions of plots (including at well locations) correspond to points with no constraining data.
564 M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
Fig. 7. Comparison of acoustic impedance estimates: (A) result of linear interpolation, (B) result of inverse distance weighting with exponent value of 1, (C) result of inverse distance weighting with exponent value of 2, and (D) result of acoustic impedance estimate obtained by Lu and McMechan (2002) by inversion of a seismic line through holes.
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 565
566
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
Fig. 8. Structural example: (A) locations of wells, vertical control lines, layer shapes and velocity data values at wells and (B) output 2D velocity model.
control horizons to define the lateral character of the structures at the required level of detail.
Acknowledgements The research leading to this paper was supported by the Sponsors of the UT-Dallas Geophysical Consortium. An initial demonstration of feasibility by Dr. Ancha Adriansyah was helpful in getting the project off the ground. Reviews by John Tipper and an anonymous reviewer improved the paper and are much appreciated. The data in Fig. 7D were kindly provided by Shaoming Lu. This paper is Contribution No. 990 from the Department of Geosciences at the University of Texas at Dallas.
References Akima, H., 1978. A method of bivariate interpolation and smooth surface fitting for irregularly distributed data points. ACM Transactions on Mathematical Software 4, 56–76. Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Elsevier, New York, 476pp. Bee, M., Jacobson, R.S., 1984. Linear inversion of body-wave data—Part III: model parameterization. Geophysics 49, 2088–2096. Berkhout, A.J., Rietveld, W.E.A., 1994. Determination of macro models for prestack migration: Part 1. Estimation of macro velocities. Proceedings of the 64th Annual International Meeting, Society of Exploration Geophysics, Expanded abstracts, Los Angeles, CA, pp. 1330–1333. Bevc, D., 1997. Flooding the topography; wave-equation datuming of level data with rugged acquisition topography. Geophysics 62, 1558–1569.
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 Bishop, T.N., Bube, K.P., Cutler, R.T., Langan, R.T., Love, P.L., Resnick, J.R., Shuey, R.T., Spindler, D.A., Wyld, H.W., 1985. Tomographic determination of velocity and depth in laterally varying media. Geophysics 50, 903–923. Cai, J., McMechan, G.A., 1999. 2-D ray-based tomography for velocity, layer shape, and attenuation from GPR data. Geophysics 64, 1579–1593. Clar, S., Ettrich, N., Ruhl, T., 1996. Can we image complex structures using smoothed macro-velocity models? Proceedings of the 66th Annual International Meeting, Society of Exploration Geophysics, Expanded Abstracts, Denver, CO, pp. 543–546. Curtis, A., Snieder, R., 1997. Reconditioning inverse problems using the genetic algorithm and revised parameterization. Geophysics 62, 1524–1525. Deutsch, C.V., Journel, A.G., 1992. GSLIB, Geostatistical Software Library and User’s Guide. Oxford University Press, New York, NY, 340pp. Dillon, W.P., Hutchinson, D.R., Drury, R.M., 1996. Seismic reflection profiles on the Blake Ridge near Sites 994, 995 and 997. In: Miller, A.T. (Ed.), Proceedings of the Ocean Drilling Program, Initial Report 164, Texas A&M University, College Station, TX, pp. 47–56. F.alt, L.M., Henriquez, A., Holden, L., Tjelmeland, H., 1991. MOHERES, a program system for simulation of reservoir architecture and properties. Proceedings of the Sixth European Improved Oil Recovery Symposium, pp. 27–39. . R., Egeland, O., 1990. SESIMIRA—a new geological Gundeso, tool for 3-D modeling of hetrogeneous reservoirs. In: Buller, A.T., Berg, E., Hjelmeland, O., Kleppe, J., Torsoeter, O., Aasen, J.O. (Eds), North Sea Oil and Gas Reservoirs—II: Proceedings of the Second North Sea Oil and Gas Reservoir Conference, Graham and Trotman, London, UK, pp. 363–371. Helgelsen, J., Magnus, I., Prosser, S., Saigal, G., Aamodet, G., Dolberg, D., Busman, S., 2000. Comparison of constrained sparse spike and stochastic inversion for porosity prediction at Kristin Field. The Leading Edge 19, 400–407. Hohn, M.E., 1999. Geostatistics and Petroleum Geology, Kluwer Academic Publishers, Dordrecht, Netherlands, 235pp. Kabir, M.M.N., Verschuur, D.J., 2000. A constrained parametric inversion for velocity analysis based on CFP technology. Geophysics 65, 1210–1222. Lu, S., McMechan, G.A., 2002. Estimation of gas hydrate and free gas saturation, concentration, and distribution from seismic data. Geophysics 67, 582–593. Mallet, J.L., 1992a. gOcad: A computer aided design program for geological applications. In: Turner, K. (Ed.), Three-
567
Dimensional Modeling with Geoscientific Information System. Kluwer Academic Publishers, Dordrecht, Netherlands, pp. 123–141. Mallet, J.L., 1992b. Discrete smooth interpolation in geometric modeling. Computer-Aided Design 24, 177–191. Mallet, J.L., 2002. Geomodeling. Oxford University Press, New York, NY, 599pp. Paull, C.K., Matsumoto, R., Wallace P.J., 1996. Leg 164 Scientific Party, Proceedings of the Ocean Drilling Program. Initial Report, 164, 1–318. Rietveld, W.E.A., Berkhout, A.J., 1994. Determination of macro models for prestack migration: Part 2. Estimation of macro boundaries. Proceedings of the 62nd Annual International Meeting, Society of Exploration Geophysics, Expanded Abstracts, Los Angeles, CA, pp. 1334–1337. Sambridge, M., Braun, J., McQueen, H., 1995. Geophysical parameterization and interpolation of irregular data using natural neighbors. Geophysical Journal International 122, 837–857. Schultz, P., Canales, L., 1997. Seismic velocity model building: CE in Dallas. The Leading Edge 16, 1063–1064. Sibson, R., 1981. A brief description of natural neighbour interpolation. In: Barnett, V. (Ed.), Interpolating Multivariate Data. Wiley, New York, pp. 21–36. Torres-Verdin, C., Victoria, M., Merletti, G., Pendrel, J., 1999. Trace-based and geostatistical inversion of 3-D seismic data for thin-sand delineation: an application in for San Jorge Basin, Argentina. The Leading Edge 18, 1070–1077. Vesnaver, A.L., 1996. Irregular grids in seismic tomography and minimum-time ray tracing. Geophysical Journal International 126, 147–165. Weber, Z., 2001. Optimizing model parameterizing in 2D linearized seismic traveltime tomography. Physics of the Earth and Planetary Interiors 124, 33–43. Wu, W.-J., Lines, L., Burton, A., Lu, H.-X., Zhu, J., Jamison, W., Bording, R.P., 1998. Prestack depth migration of an Alberta Foothills data set—the Husky experience. Geophysics 63, 392–398. Xu, W., Tran, T.T., Srivastava, E.M., Journel, A.G. 1992. Integrating seismic data in reservoir modeling: the colocated cokringing alternative. Proceedings of the 67th Annual Technical Conference of the Society of Petroleum Engineers, Washington, Article SPE 24742, pp. 833–842. Zhu, J., Lines, L.R., 1998. Comparison of Kirchhoff and reverse-time migration methods with applications to prestack depth imagings of complex structures. Geophysics 63, 1166–1176.
Building 2-D stratigraphic and structure models from well log data and control horizons Maulin D. Patel, George A. McMechan* Center for Lithospheric Studies, University of Texas at Dallas, P.O. Box 830688, 2601 N. Floyd Road, Richardson, TX 75083-0688, USA Received 25 April 2002; received in revised form 3 December 2002; accepted 15 December 2002
Abstract An algorithm to build a gridded 2-D seismic velocity (or any other physical property) model from well log data and control horizons is developed. Interpolation of well log data onto a 2-D grid uses inverse distance weighting or linear interpolation, guided by the shape of the control horizons that are predefined from seismic or other 2-D constraints. A key feature of the models is that they may contain layers that are truncated by unconformities or at faults, or that lap out smoothly at their tops and bottoms. Abrupt or smooth terminations are controlled by user flags. Applications are illustrated using resistivity and acoustic impedance log data from the Blake Ridge (offshore Carolina) and using a complex structure produced by overthrust tectonics in western Canada. Geologically reasonable models can be produced only if there are sufficient wells to sample every salient element in the model and sufficient control horizons to define the lateral character of the structures at the required level of detail. r 2003 Elsevier Science Ltd. All rights reserved. Keywords: Model building; Stratigraphic interpolation; Well log interpolation; Structural control
1. Introduction Time-to-depth conversion, numerical modeling, and depth imaging of seismic data require geologically and geophysically realistic velocity models (Berkhout and Rietveld, 1994; Rietveld and Berkhout, 1994). Increasing quality of seismic data and survey accuracy demand correspondingly high-quality velocity models. Models can be parameterized and interpolated in a variety of ways. Tomographic models are usually either constant property layers (Cai and McMechan, 1999) or regular grids (Bishop et al., 1985); recent advances include irregular grids that adapt with finer sampling as data resolution increases (Sambridge et al., 1995; Vesnaver, 1996; Curtis and Snieder, 1997; Weber, 2001). Withinlayer gradients are also possible (Kabir and Verschuur, 2000). A number of sophisticated commercial 3-D model builders also exist (Schultz and Canales, 1997; Mallet, *Corresponding author. Tel.: +972-883-2424; fax: +1-972883-2829. E-mail address: [email protected] (G.A. McMechan).
1992a, 2002; F.alt et al., 1991). Most of these are directed toward construction of either macro models for input to imaging (Clar et al., 1996), or detailed petrophysical models as part of seismic inversion (Bee and Jacobson, 1984), or for input to reservoir simulations (Aziz and Settari, 1979; Gundeso. and Egeland, 1990). Other applications include modeling stratigraphic channels and detecting faults (Mallet, 2002). Model building explicitly requires interpolation to provide consistent parameter values between locations where there are measurements. In the past, this has usually been done with statistical models and various forms of kriging (e.g. Hohn, 1999; Xu et al., 1992; Deutsch and Journel, 1992; Mallet, 2002). However, where deterministic constraints are available (for example, from seismic or ground penetrating radar data), it is possible to produce models with deterministic geometry. Examples of comparisons of statistical and deterministic models are given by Torres-Verdin et al. (1999) and Helgelsen et al. (2000). Where only widely spaced well logs are available, it is impossible to adequately characterize lateral changes in lithology because of
0098-3004/03/$ - see front matter r 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0098-3004(03)00039-6
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
558
insufficient sampling. This limits the features that can be reliably interpolated to those of the largest scales. Where fairly continuous profiles are available (such as seismic data) the observed structural geometry can be inserted as constraints during interpolation of the well log data. Many deterministic interpolations (including for irregular grids) already exist (e.g., Akima, 1978; Sibson, 1981; Mallet, 1992b). These tend to be rather complicated, and the commercial versions are very expensive. Thus, it seems worth to extract only the essential concepts and implement them in a small, flexible, easyto-use format. This paper describes a procedure for building 2-D velocity (or other physical property) models, given well log data and control horizon geometries. The control horizons are typically obtained from migrated seismic depth sections. Measurements of some physical property (say seismic velocity, density, resistivity or porosity) constitute the well log data. Well data provide depths to key formation boundaries from drilling data or from well logs. The algorithm interpolates well log data so as to find reasonable values of these properties at all points in a 2-D grid by following the spatial geometry of the control horizons. Linear interpolation or inverse distance weighting are used for interpolation. Linear interpolation uses data only from the wells to either side of a specified location; inverse distance weighting involves contributions from all wells. A variety of layer truncation types are included, and the algorithms can be extended to 3-D.
2. Model description A vertical earth section is represented on a 2-D plane (Fig. 1). The control horizons are lines that go from one side of the model to the other to define the layer shapes; these will generally be a combination of sedimentologi-
X
cal bounding surfaces and other internal features. Vertical lines are the positions of wells that provide the well log data to be interpolated. Interpolation is between the well logs, using ties between the wells specified by the predefined control horizons and flags that specify how the interpolation is to be done. The final output is a 2-D grid of interpolated values. 2.1. Model specifications The number and positions of the vertical model control lines are arbitrary, but they should be placed to provide accurate descriptions of the model geometry. The X positions of the vertical control lines include the locations of all the wells. The algorithm has an option of omitting any well(s) from the model interpolation. The model can extend beyond the leftmost and rightmost X positions of the vertical control lines. The depth of each control horizon is specified on each vertical control line. The uppermost control horizon (representing the air– earth interface) can be flat or curved and is not necessarily at depth ¼ 0: The input well data points start from control horizon # 1 (not necessarily from depth ¼ 0). The depth extents of the well data are arbitrary. The output grid spacing in the vertical direction (DZ) is equal to the vertical spacing between samples in the well log data, and is the same in all well logs. All control horizons extend from left to right across the entire section. Truncated model features are defined by curves that also cross the entire section, but coincide with other control horizons beyond points where the model features are truncated. The control horizons are numbered from top to bottom in increasing order (Fig. 1). Control horizons are single-valued in depth and do not cross each other. Effectively multivalued geological horizons are composed by connecting
Position Well #1
Well #2
Well #3
Depth
Z
Control Horizon 4
Control Horizons
Control Horizon 5
Fig. 1. A vertical earth section on a 2-D plane with horizontal coordinate X and vertical coordinate Z: Control horizons and well logs provide geological constraints and geophysical data to be interpolated, respectively.
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
segments of single-valued horizons. The geometry of the control horizons conform to the geological environment and may, for example, be provided by seismic reflection profiles.
3. The algorithm The parameters input to the algorithm are: the number of output grid points in the X (horizontal) and Z (vertical) direction, the maximum number of samples in the well log data among all the wells at a depth increment equal to the vertical grid spacing (DZ), the grid spacing in the horizontal direction (DX ), the grid spacing in the vertical direction (DZ), the number and positions of the well(s) in the X direction, the amount of smoothing to apply to the well logs, the number and coordinates of vertical control lines, the boundary treatment parameters for all the control horizons (described below), and the well log data for all the wells. The task is to interpolate the well log data to generate consistent values of the measured properties at all points in a 2-D grid by following the geometry defined by the control horizons. Linear interpolation or inverse distance weighting can be chosen for the interpolation. The output is a 2-D gridded distribution of the estimated physical parameter on a vertical section through the earth. Step-by-step execution of the algorithm is explained below. 3.1. Outline of the algorithm The first step is to compute the depths of the control horizons at each of the vertical control lines, which include (but are not limited to) all of the wells. When the X coordinate of a control point coincides with a well position, the Z coordinates constitute the control horizon depths at that well. Where a geological feature has pinched out or has been truncated, the corresponding control horizons are superimposed on other control horizons that do extend completely across the model; thus, two or more control horizons may be at the same depth at a well, but are not allowed to cross. In Fig. 1, wells #1 and #3 do not have grid points between control horizons 4 and 5, so if well #2 is omitted from consideration, interpolation cannot be performed between horizons 4 and 5. If the number of grid points between a pair of control horizons at a particular well is zero, then the corresponding well data are assigned the value 1: When the control horizon #1 is not at depth ¼ 0; grid points are inserted above first control horizon. These points will not have any data values from wells and so are also assigned the value of 1: Similarly, when a well does not extend to the grid bottom, data points inserted below the bottom of the well are also assigned the value 1: Thus, the processed well data cover the
559
entire depth of the model, whereas the actual well data may not. When a data value of 1 is found during the interpolation, that data value is ignored as if that well was not present. The algorithm works from top to bottom and from left to right through all the grid points. For each horizontal position of the grid point, there are four possibilities: the grid point is to the left of the leftmost well position, or the grid point is to the right of the rightmost well position, or it superimposes on a well position, or it is between two well positions. For inverse distance weighting the first, second and fourth possibilities are treated in the same way and an inverse distance weighting routine is called. Where a grid point superimposes on a well position, the data values used are those of the well data at that location (no interpolation). When linear interpolation is specified, a linear interpolation routine is called. All the grid points to the left of the leftmost well are assigned the data values of that well (no interpolation). Similarly all the grid points to the right of the rightmost well are assigned the well data values at the same stratigraphic level in that well (no interpolation). For the third and fourth possibilities, the linear interpolation routine finds the interpolated data values. To interpolate a value at point a (Fig. 2A), the values at point b in well #1 and point c in well #2 are used; a–c are all at the same relative vertical position between control horizons j 1 and j: For example, if a is 0.3 of the increment from j 1 to j; then we can find the data values at b and c; which are also 0.3 of the increment from j 1 to j; at wells #1 and #2. This is done by subsampling the well data so there are the same number of points between the two given horizons at all wells. Layer truncations are a special case that is described in Section 3.4. 3.2. Linear interpolation Consider linear interpolation for the representative point a in Fig. 2A. The value assigned at point a is given by Va ¼
Da Db ðVc Vb Þ þ Vb ; Dc Db
where Va ; Vb and Vc are the data values at points a; b; and c; respectively. Da ; Db ; and Dc are the distances of points a; b and c from the origin, respectively. 3.3. Inverse distance weighted interpolation Consider the representative point d on the dashed interpolation trajectory in Fig. 2B. The value assigned to the point d is given by P Wi Vi Vd ¼ Pi ; i Wi
560
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 Position Control horizon j-1
Depth
Well #1
c a
b
Control horizon j
(A) Well #3
Depth
Well #2
Well #4
Well #5
C ontrol Region 1 e
g Region 2
f
d
(B) Fig. 2. Geometry for (A) linear interpolation and (B) weighted interpolation. In (A) dashed line maintains same relative position between control horizons j 1 and j: In (B) wells #3 and #5 have no data values for region 1, so for all grid points in region 1, only data values of well #4 are considered in interpolation. Interpolation at point d; on dashed interpolation trajectory efg, uses data values at f and g for linear interpolation, or at e–g for inverse weighted interpolation. Interpolations of grid points to left of well #3 and to right of well #5 use data values from only well #3 and well #5, respectively, for linear interpolation, or weighted values from all wells, for inverse weighted interpolation.
where the sum is over the contributions from each of the available wells (e.g. i ¼ 3; 4; 5). Vi is the data value in the well i at the same relative position that point d is between the control horizons above and below d; and Wi is the weight given to Vi : The sum in the denominator normalizes the weights so that their sum is 1.0: Wi ¼ ðDid Þp ; where Did is the distance from well i to point d; and p is the user-selected power weighting of the contributions. To avoid division by zero, the observed data values are used at the wells (where Did ¼ 0). A data value of 1 is ignored as if that well was not present in the 2-D grid, so the weight of that well is not calculated and it will have no impact on the weighted value. 3.4. Truncated control horizons All the control horizons extend from left to right across the entire section. To handle the problem of smaller curve segments, which do not cross the entire section, they are superimposed on control horizons that do extend from left to right. For example, in Fig. 3A, horizon db ends at b; but it needs to be specified as dbgc. Similarly, horizon eg ends at g; but is specified as efgc. With the interpolation described above, if a control horizon is truncated, then the interpolated trajectories will converge at a point. For example, in Fig. 3A, the dashed interpolated trajectories between the input
control horizons dbgc and efgc will all converge at point g: If bg is a stratigraphic pinchout that occurred during the initial deposition of the layer, then this convergent interpolation geometry is geologically reasonable. However, it is also necessary to consider the alternative that horizon abc is an erosional unconformity that has truncated horizons db and eg; in this case, the appropriate geometry is obtained by extrapolation past b; guided by the single existing horizon fg (Fig. 3B). The data values along the vertical line segment bf are replicated vertically at each horizontal position along the curve fg (see the dashed extrapolation trajectories in Fig. 3B; values extrapolated above bc are discarded). Layers that are truncated at their bottom end, for example by a fault, can be similarly extrapolated to produce the desired corresponding geometry. Thus, there are three possible behaviors that the user can specify: layer convergence, layer truncation at its upper end, or layer truncation at its bottom end. The decision on which choice is appropriate should come from a priori information from the geological environment. A similar, but more complicated solution to the layer truncation problem is illustrated by Mallet (2002); this involves one real and one virtual control horizon.
4. Examples In this section we present a series of examples from both synthetic and field data that illustrate the
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 Position
561
Position h
Depth
a
c
b
A
c
b g
g d
D f
f e
(A)
E
(B)
Fig. 3. Handling of truncated control horizons. (A) With layer convergence, dashed interpolated trajectories between input control horizons dbgc and efgc will converge at point g: (B) With truncation, appropriate geometry is obtained by extrapolation, guided by single existing horizon fg.
capabilities of the software. The examples show the behavior of many of the options described above. We show interpolation of resistivity and acoustic impedance well logs from the Blake Ridge, and compare with the results from a seismic inversion. This region is stratigraphically complicated, but structurally simple. Then, we apply the software to building a model from a structurally complex overthrust region in western Canada.
4.1. Synthetic example This example uses synthetic data and control horizons that illustrate the range of capabilities of the software. Thus, this model is not constrained to being geologically reasonable. Fig. 4A shows the locations of 16 vertical control lines (three of which are wells), 24 control horizons, and the range of data values at the wells. The boundary treatment parameters for control horizons and the interpolation algorithms to be used are specified by the user. Fig. 4B shows the estimates obtained by linear interpolation. Note the ability to deal with layer convergence (at a), truncation of a layer at its upper end (at b), and layer truncation at its bottom end (at c). Weighted interpolation (Fig. 4C) provides estimates when the X positions of the control points do not extend across the entire model. Removing the first and last vertical control lines results in the control horizons being extended horizontally to the model edges using the same Z positions as at the leftmost and rightmost input vertical control lines, respectively. Fig. 4D shows the effect of changing boundary treatment parameters from truncation to smooth convergence everywhere, and omitting well #2 from consideration. Omission of well #2 results in region d (which gets data only from well #2) turning white (i.e. becomes undefined), and specifying layer convergence of all boundaries results in all interpolation trajectories converging at points. Thus, knowing either layer shapes or having logs of the values to be interpolated is not sufficient; both must be present.
The interpolation procedures described above cannot by themselves produce geologically valid models. The models are only as complete as the constraints that are provided. For example, if only two well logs that are 5 km apart are given, the interpolation necessarily produces continuity across B2:5 km (if the quantity being interpolated has different values at the corresponding relative depths in the wells) or across all 5 km (if the quantity log interpolated has similar values at the corresponding relative depths). Additional lateral detail can be inserted only by providing additional constraints, for example, on the geometry (perhaps from seismic control) between the holes (e.g. Fig. 4) or by requiring the interpolation to confirm to a statistical model that defines a distribution of expected lateral dimensions. Thus, this approach has the advantage that the result is faithful to the type and resolution of the information that is input.
4.2. Stratigraphic examples The examples in this section are based on well logs from Ocean Drilling Program (ODP) Leg 164 holes 994D, 995B, and 997B on the Blake Ridge off the east coast of North America. Details of the data acquisition and processing are presented by Paull et al. (1996), Dillon et al. (1996), and Lu and McMechan (2002). Fig. 5 shows the locations of the three holes, the five control horizons that were extracted (as coherent reflections) from a seismic line through the holes, and the resistivity well log profiles measured in the holes. With these as input, we construct the interpolated sections in Fig. 6 using four different interpolations. The main differences are near the edges of the section. Linear interpolation (Fig. 6A) includes extrapolation, before the first and beyond the last holes, with the same data values as at those holes and so provides estimates only over the stratigraphic depths at which the logs exist (compare with Fig. 5). Weighted interpolations (Figs. 6B–D) provide estimates over all stratigraphic
Fig. 4. Synthetic example: (A) Locations of wells, control horizons, vertical control lines and data values at wells. Numbers x þ y þ z represent velocities at that well between a pair of control horizons in order of increasing depth. (B) Result of linear interpolation. (C) Result of inverse distance weighting with exponent value of 2 and omitting first and last vertical control lines. (D) Result of inverse distance weighting with exponent value of 2 and omitting well #2; specifying convergent layer boundaries results in layer convergence everywhere.
562 M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
563
Fig. 5. Stratigraphic example. Figure shows locations of wells, control horizons and resistivity log profiles measured in wells of ODP Leg 164 holes 994D, 995B, and 997B on Blake Ridge off east coast of North America. Corresponding interpolated sections are in Fig. 6.
depths at which a log exists in any hole. Away from the holes, the weighted extrapolations approach the average value of the logs, at the corresponding stratigraphic depth, at a rate that depends on the chosen exponent; for example, at the right side of Fig. 6D (with exponent 2), the rate of change is more rapid that in the same region of Fig. 6B (with exponent 0). Figs. 7A–C show the estimates obtained by linear and weighted interpolations of the acoustic impedance (the product of P-wave velocity and density) logs using the same control horizons as used in Fig. 6. For comparison, Fig. 7D shows the acoustic impedance estimated by inversion of the seismic section between the holes (after Lu and McMechan, 2002). The overall behaviors of the various estimates are similar; in particular, the low-impedance layer near 550–600 m depth (which also has high resistivity in Fig. 6) is clearly visible in both the estimates obtained by interpolation and by seismic inversion. This layer was interpreted by Lu and McMechan (2002) as containing free methane. 4.3. Structural example In addition to stratigraphic interpolation/extrapolation, as illustrated in the foregoing subsection, the software has utility in building structurally complicated (multivalued) models. Fig. 8A contains the vertical control positions and the layer shapes for building a velocity model for a data set used in SEG workshops on structural imaging in 1995 and 1996. Results based on this model have been presented by many authors, including Bevc (1997), Wu et al. (1998), and Zhu and Lines (1998). Fig. 8B contains the output 2-D velocity model. Note the ability to deal with multivalued
horizons (e.g., in region a of Fig. 8B), including faults (e.g., in region b of Fig. 8B), layer truncations (e.g., in region c of Fig. 8B), and pinchouts (e.g., in region d of Fig. 8B). Thus, all the main types of structural relationships can be included. Although Fig. 8A appears complicated, its generation is facilitated by a graphical user interface that allows points to be added, moved or deleted on a computer monitor, with mouse commands. This makes use in subsequent applications (e.g., updating velocities or layer geometries in iterative prestack migration) not only feasible, but also user-friendly. Extension to 3-D model building is also feasible, requiring mainly that the control horizons be specified as 3-D surfaces rather than as 2-D lines.
5. Summary We have developed and illustrated the operation of an algorithm for building 2-D models by interpolation and extrapolation of well log data. The interpolations and extrapolations are constrained by control horizons that are defined by seismic (or other) data that provide geometrical guides between the holes. Interpolation can be linear (between holes) or inverse distance weighting from all holes. Examples include interpolation between logs from a stratigraphically complex, but structurally simple, region off the east coast of North America, and building a structurally complex model for an overthrust zone in eastern Canada. Geologically reasonable models can be produced only if there are sufficient wells to sample every salient element in the model and sufficient
Fig. 6. Comparison of resistivity estimates: (A) result of linear interpolation, (B) result of inverse distance weighting with exponent value of 0, (C) result of inverse distance weighting with exponent value of 1, and (D) result of inverse distance weighting with exponent value of 2. White portions of plots (including at well locations) correspond to points with no constraining data.
564 M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
Fig. 7. Comparison of acoustic impedance estimates: (A) result of linear interpolation, (B) result of inverse distance weighting with exponent value of 1, (C) result of inverse distance weighting with exponent value of 2, and (D) result of acoustic impedance estimate obtained by Lu and McMechan (2002) by inversion of a seismic line through holes.
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 565
566
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567
Fig. 8. Structural example: (A) locations of wells, vertical control lines, layer shapes and velocity data values at wells and (B) output 2D velocity model.
control horizons to define the lateral character of the structures at the required level of detail.
Acknowledgements The research leading to this paper was supported by the Sponsors of the UT-Dallas Geophysical Consortium. An initial demonstration of feasibility by Dr. Ancha Adriansyah was helpful in getting the project off the ground. Reviews by John Tipper and an anonymous reviewer improved the paper and are much appreciated. The data in Fig. 7D were kindly provided by Shaoming Lu. This paper is Contribution No. 990 from the Department of Geosciences at the University of Texas at Dallas.
References Akima, H., 1978. A method of bivariate interpolation and smooth surface fitting for irregularly distributed data points. ACM Transactions on Mathematical Software 4, 56–76. Aziz, K., Settari, A., 1979. Petroleum Reservoir Simulation. Elsevier, New York, 476pp. Bee, M., Jacobson, R.S., 1984. Linear inversion of body-wave data—Part III: model parameterization. Geophysics 49, 2088–2096. Berkhout, A.J., Rietveld, W.E.A., 1994. Determination of macro models for prestack migration: Part 1. Estimation of macro velocities. Proceedings of the 64th Annual International Meeting, Society of Exploration Geophysics, Expanded abstracts, Los Angeles, CA, pp. 1330–1333. Bevc, D., 1997. Flooding the topography; wave-equation datuming of level data with rugged acquisition topography. Geophysics 62, 1558–1569.
M.D. Patel, G.A. McMechan / Computers & Geosciences 29 (2003) 557–567 Bishop, T.N., Bube, K.P., Cutler, R.T., Langan, R.T., Love, P.L., Resnick, J.R., Shuey, R.T., Spindler, D.A., Wyld, H.W., 1985. Tomographic determination of velocity and depth in laterally varying media. Geophysics 50, 903–923. Cai, J., McMechan, G.A., 1999. 2-D ray-based tomography for velocity, layer shape, and attenuation from GPR data. Geophysics 64, 1579–1593. Clar, S., Ettrich, N., Ruhl, T., 1996. Can we image complex structures using smoothed macro-velocity models? Proceedings of the 66th Annual International Meeting, Society of Exploration Geophysics, Expanded Abstracts, Denver, CO, pp. 543–546. Curtis, A., Snieder, R., 1997. Reconditioning inverse problems using the genetic algorithm and revised parameterization. Geophysics 62, 1524–1525. Deutsch, C.V., Journel, A.G., 1992. GSLIB, Geostatistical Software Library and User’s Guide. Oxford University Press, New York, NY, 340pp. Dillon, W.P., Hutchinson, D.R., Drury, R.M., 1996. Seismic reflection profiles on the Blake Ridge near Sites 994, 995 and 997. In: Miller, A.T. (Ed.), Proceedings of the Ocean Drilling Program, Initial Report 164, Texas A&M University, College Station, TX, pp. 47–56. F.alt, L.M., Henriquez, A., Holden, L., Tjelmeland, H., 1991. MOHERES, a program system for simulation of reservoir architecture and properties. Proceedings of the Sixth European Improved Oil Recovery Symposium, pp. 27–39. . R., Egeland, O., 1990. SESIMIRA—a new geological Gundeso, tool for 3-D modeling of hetrogeneous reservoirs. In: Buller, A.T., Berg, E., Hjelmeland, O., Kleppe, J., Torsoeter, O., Aasen, J.O. (Eds), North Sea Oil and Gas Reservoirs—II: Proceedings of the Second North Sea Oil and Gas Reservoir Conference, Graham and Trotman, London, UK, pp. 363–371. Helgelsen, J., Magnus, I., Prosser, S., Saigal, G., Aamodet, G., Dolberg, D., Busman, S., 2000. Comparison of constrained sparse spike and stochastic inversion for porosity prediction at Kristin Field. The Leading Edge 19, 400–407. Hohn, M.E., 1999. Geostatistics and Petroleum Geology, Kluwer Academic Publishers, Dordrecht, Netherlands, 235pp. Kabir, M.M.N., Verschuur, D.J., 2000. A constrained parametric inversion for velocity analysis based on CFP technology. Geophysics 65, 1210–1222. Lu, S., McMechan, G.A., 2002. Estimation of gas hydrate and free gas saturation, concentration, and distribution from seismic data. Geophysics 67, 582–593. Mallet, J.L., 1992a. gOcad: A computer aided design program for geological applications. In: Turner, K. (Ed.), Three-
567
Dimensional Modeling with Geoscientific Information System. Kluwer Academic Publishers, Dordrecht, Netherlands, pp. 123–141. Mallet, J.L., 1992b. Discrete smooth interpolation in geometric modeling. Computer-Aided Design 24, 177–191. Mallet, J.L., 2002. Geomodeling. Oxford University Press, New York, NY, 599pp. Paull, C.K., Matsumoto, R., Wallace P.J., 1996. Leg 164 Scientific Party, Proceedings of the Ocean Drilling Program. Initial Report, 164, 1–318. Rietveld, W.E.A., Berkhout, A.J., 1994. Determination of macro models for prestack migration: Part 2. Estimation of macro boundaries. Proceedings of the 62nd Annual International Meeting, Society of Exploration Geophysics, Expanded Abstracts, Los Angeles, CA, pp. 1334–1337. Sambridge, M., Braun, J., McQueen, H., 1995. Geophysical parameterization and interpolation of irregular data using natural neighbors. Geophysical Journal International 122, 837–857. Schultz, P., Canales, L., 1997. Seismic velocity model building: CE in Dallas. The Leading Edge 16, 1063–1064. Sibson, R., 1981. A brief description of natural neighbour interpolation. In: Barnett, V. (Ed.), Interpolating Multivariate Data. Wiley, New York, pp. 21–36. Torres-Verdin, C., Victoria, M., Merletti, G., Pendrel, J., 1999. Trace-based and geostatistical inversion of 3-D seismic data for thin-sand delineation: an application in for San Jorge Basin, Argentina. The Leading Edge 18, 1070–1077. Vesnaver, A.L., 1996. Irregular grids in seismic tomography and minimum-time ray tracing. Geophysical Journal International 126, 147–165. Weber, Z., 2001. Optimizing model parameterizing in 2D linearized seismic traveltime tomography. Physics of the Earth and Planetary Interiors 124, 33–43. Wu, W.-J., Lines, L., Burton, A., Lu, H.-X., Zhu, J., Jamison, W., Bording, R.P., 1998. Prestack depth migration of an Alberta Foothills data set—the Husky experience. Geophysics 63, 392–398. Xu, W., Tran, T.T., Srivastava, E.M., Journel, A.G. 1992. Integrating seismic data in reservoir modeling: the colocated cokringing alternative. Proceedings of the 67th Annual Technical Conference of the Society of Petroleum Engineers, Washington, Article SPE 24742, pp. 833–842. Zhu, J., Lines, L.R., 1998. Comparison of Kirchhoff and reverse-time migration methods with applications to prestack depth imagings of complex structures. Geophysics 63, 1166–1176.