ABSTRACT

Dynamic forward modeling of carbonate reservoir sequence stratigraphy and diagenetic overprint can yield rapid, cost-effective reservoir characterization. The common practice in reservoir characterization now relies heavily on massive data accumulation and geostatistics to produce the three-dimensional geocellular static model which is the basis for flow simulation. In dynamic forward modeling, reliance on understanding of geological process allows high resolution prediction of the geometry of permeable and impermeable units and horizons within the reservoir. Data requirements are reduced to state-of-the-art information on a relatively small number of control wells which constrain and calibrate the forward model. Sensitivity-testing among formally-stated competing concepts is encouraged. In the long-term, it is the accurate prediction of reservoir response to future production that will afford choice among competing static models and flow simulations. The goal should be to predict future problems and avoid them, rather than wait to observe problems and react to them.

INTRODUCTION

The advent of user-friendly computer work-stations continues to bring change to the daily activities of petroleum geoscientists and reservoir engineers. Thus far, the new technology has been used to do more of the same, but to do it more efficiently. Computer memory holds data that would previously have occupied several rooms of filing cabinets. Data retrieval time is orders of magnitude more rapid. Geologists are correlating more and more electric logs, describing more and more core, and adding more names to the list of rock-types they track. Likewise, a flow calculation that used to take weeks with a calculator is now accomplished in minutes on a computer. There remains the scarcely-tapped opportunity to engage the computer in dynamic forward modeling toward a whole new level of rapid, cost-effective carbonate reservoir characterization and flow simulation.

In dynamic forward modeling, the paradigm or “framework within which we attempt to solve problems” is separable into topics that can be studied independently. For carbonate reservoir geology, these separable topics include at least the following frameworks: tectonic/isostatic, eustatic (including both tectono-eustatic and glacio-eustatic components), climatic, sedimentologic, and diagenetic. Instead of attempting to unravel observed “detail complexity”, as is the tradition of the carbonate reservoir geologist, “dynamic complexity” is created by combining processes within the computer-based dynamic forward model. The output of the forward model is then compared to reality. When the complexity created resembles the complexity observed, a new understanding of the carbonate reservoir is at hand.

This paper reviews the promise and the pitfalls of computer-based dynamic forward modeling with particular reference to Arabian Gulf Mesozoic carbonate reservoir sequence stratigraphy and diagenetic overprint. The underlying philosophy of dynamic forward modeling is reviewed with particular reference to current and future practices in carbonate reservoir characterization, flow simulation and reservoir-performance forecasting (e.g. Saleri, 1998).

MODEL CONSTRUCT

Numerous, relatively simple processes can be brought together in a dynamic forward model to produce the seeming complexity which is observed in carbonate reservoir studies. Taken individually or in pairs, the relationships involved are easy to understand. Tectonic subsidence and sea-level fluctuation combine to produce accommodation space for new sediments or subaerial exposure of previously deposited sediments. The mud content of carbonate sediments is often a function of water depth. Fresh-water diagenesis is quite different above and below the top of the water table, etc. Yet when nature presents all of this at once in a sequence of rocks, the complexity is more than the human mind can easily comprehend. The dynamic forward model becomes an extension of the mind, allowing the problem to be dealt with efficiently.

Figure 1 presents a flow diagram which progresses from simple generalities above to complex combinations below. In a forward model, these relationships are depicted as equations. While very few things in geology are truly independent, the math can specify that independent variables change only as a function of time whereas dependent variables change as a function of two or more independent variables. For example, subsidence and sea-level are here considered independent variables whereas accommodation space changes as the combination of the two.

Most of the topics discussed below are thriving fields of science unto themselves. The science is available for the forward model to be made as complex as may be required by the task at hand. However; the general rule recommended here is to construct the simplest, geologically reasonable model that is consistent with the target data.

Tectonics

Basin subsidence is commonly the most important independent variable in the generation of regional stratigraphy. Why and when did the basin subside? Was the process continuous or episodic? Over what scale was the load isostatically compensated? The answers to these questions also serve to constrain geological reasonableness with regards to trade-offs among subsidence history, sea-level history, and geochronologic uncertainty.

Plate Tectonics Context

With regards to vertical motion of the lithosphere, subsidence driven by lithospheric cooling is by far the best understood process. The fundamental observations come from the study of oceanic lithosphere (e.g. Parsons and Sclater, 1977; Renkin and Sclater, 1988), but the physics can be applied generally. The basic observation is that the depth to the top of oceanic lithosphere varies as a function of the age-of-oceanic lithosphere. The hot, newly-formed oceanic lithosphere starts off about 2.8 kilometers (km) below sea-level. Cold, old oceanic lithosphere resides beneath 5 to 6 km of water. As the lithosphere cools, it becomes more dense and the accompanying loss of volume produces this relationship.

This relationship extends very nicely to continental passive margins. Continental lithosphere is reheated prior to continental break-up. Subsidence ensues following continental break-up and formation of passive margins around a newly-opening ocean basin (e.g. Watts and Thorne, 1984; Lister et al., 1986; Braun and Beaumont, 1989). Subsidence rate varies as a function of lithosphere thickness (Fowler, 1990, p. 378-398). Figure 2 presents a typical curve for subsidence rate as a function of time since passive margin formation. The passive margin model is probably sufficient for much of Arabian Gulf Mesozoic carbonate sedimentation. The predicted subsidence history is non-linear.

Subsidence and Isostasy

The aerial extent of isostatic compensation in response to new load, is to a first approximation, a function of the flexural rigidity of the lithosphere, which is, in turn, controlled primarily by its thickness (Forsyth, 1985). Note well that variation in the flexural rigidity of the lithosphere may be expected across the length of a geologic cross-section. For example, passive margins may start out in continental lithosphere, go next to fractured continental lithosphere, and finally to oceanic lithosphere. Each of these will have different wave length and relaxation time parameters (e.g. Turcotte and Schubert, 1982).

Interestingly, whereas the new weight of sediment is isostatically compensated substantially on the time scale of high-stands and low-stands in the model, diagenetic alteration associated with phreatic lenses can almost certainly be regarded as requiring no further compensation. The calcium carbonate is already there as sediment; diagenesis mostly just redistributes it. Thus, the thickness of higher-order high-stand sequence units reflects sea-level amplitude plus isostasy. The thickness between top of third-order high-stand units and top of third-order low-stand phreatic lenses reflects only sea-level amplitude.

Eustasy

There are pitfalls with regards to sea-level estimates derived internally from stratigraphic sequences under consideration. Where possible, independent estimates of the tectonic and ice volume components which affect sea-level time series should be utilized. For time intervals where such data are unavailable, sea-level may be inferred from similar intervals on the basis of process/response models.

Internally-derived Sea-level Estimates

It has long been popular to attempt to derive sea-level history directly from the stratigraphic sequences themselves. In part, this was driven by necessity, in that interest in sea-level history preceded tectono-eustatic and glacio-eustatic technology. In part, this was simply the traditional inductive mode, “Let the data speak for itself”, which has driven stratigraphic thinking for years. Viewed from the perspective of process, the geological reasonableness of internally-derived sea-level estimates is questionable.

“Fischer plots” (Fischer, 1964) remain a popular internally-derived sea-level estimate. The Fischer construct assumes a long-term trend of subsidence and/or sea-level (i.e., local relative sea-level) based upon environmental datums at bottom and top of the section. Then, internal stratigraphy is subdivided by assuming that each cycle represents a constant period. Thus varying thicknesses of cycles is taken to represent greater or lesser heights of the various stands of the sea. The Fischer plot is a numerologically satisfactory description of the sedimentary cycles observed in the stratigraphy. Unfortunately, this construct has no claim to geological reasonableness, because it relies on assumptions which are too simple regarding periodicity of cycles.

Modern sequence stratigraphy (e.g. Vail et al., 1977; Vail and Hardenbol, 1979; Haq et al., 1987; Moore et al., 1987; Greenlee and Moore, 1988; Posamentier et al., 1988) builds upon the old physical stratigraphy concept of continental freeboard (originally a nautical term; with no freeboard, a ship is submerged) (e.g. Sloss, 1963, 1964). The basis for the modern freeboard construct is the regional seismic line and/or cross-section constructed perpendicular to strike and away from local tectonic complications. A series of relative sea-level changes is proposed to explain the geometric relationships within the regional stratigraphy. After compiling these reconstructions for the same time interval from several regions, the similarities are taken as a global eustatic signal. Such a construct for the last 65 million years (My) is presented in the left portion of Figure 3. Later in this paper, this construct is compared to independent estimates of tectono-eustasy and glacio-eustasy. For now, suffice it to say the curve is of low resolution and depicts sea-level amplitude about three times greater than any currently recognized geologically reasonable mechanism would produce.

Both the Fischer plot and the sequence stratigraphy generalized curve camouflage significant pitfalls on the path to quantitative stratigraphic problem-solving. Both are commonly taken to represent estimates of sea-level fluctuation, when both are in fact geologically unreasonable; the Fischer plot because of its assumption of equal cycle duration; the sequence stratigraphy generalized curve because of its low resolution and excess amplitude. Both may be qualitatively useful in hands-on geology; but, as forward modeling sea-level input files, geological reasonableness is required.

Independent Estimates of Tectono-eustatic Sea-level History

If the volume of the global ocean basin changes and the water volume and area of the Earth’s surface remain constant, then sea-level relative to the continents must rise or fall. Pitman (1978), Kominz (1984), and Harrison (1988) evaluate the changing volume of the global ocean basin over the last 80 My. The volume changes as a function of sea-floor spreading rate, subsea volcanism, changing aerial extent of continents caused by continental collision/stretching, and sediment accumulation in the marine environment (Table 1). By quantitatively evaluating the importance of each process, the effect on continental freeboard can be evaluated (Figure 4). Values and resultant curves are undoubtedly subject to further refinement. For this discussion, the point is that each of these processes can be studied independently and their effects upon sea-level quantified through time. It should not be necessary to rediscover the tectono-eustatic effect in each and every new local study of a carbonate sequence.

Several important conclusions with regards to the entire Phanerozoic can be drawn from the data presented in Table 1 and Figure 4. First and foremost, these effects produce sea-level changes that occur relatively slowly. Inspection of Table 1 indicates 5 millimeters per thousand years (mm/Ky) is an extremely rapid rate of change for these tectono-eustatic mechanisms. Note also that the 80 My curve of Harrison (1988) bears strong resemblance to the corresponding portion of Sloss (1963) generalizations regarding continental freeboard throughout the Phanerozoic. While continental freeboard varies through time from continent to continent (Bond, 1979, 1985; Harrison et al., 1983, 1985; Harrison, 1988), it is reasonable that processes and rates from the study of the last 80 My can be extrapolated throughout the Phanerozoic.

The Oxygen Isotope Record of Glacio-eustasy

The oxygen isotope composition of deep sea planktic and benthic foraminifers includes information regarding glacio-eustatic sea-level history. These studies of relatively young sediments serve as calibration data for orbital forcing of glacio-eustasy calculations review which follows.

As ice accumulates on the continents, the ice stores isotopically light water, making the well-mixed world ocean isotopically heavier with respect to Oxygen-18. A sizable literature indicates the Late Pleistocene oxygen isotope record is a good approximation of sea-level history (e.g. Fairbanks and Matthews, 1978; Bender et al., 1979; Matthews, 1986, 1990). The amplitude of Late Pleistocene glacio-eustatic sea-level fluctuation is on the order of 100 meters (m). This implies rates of sea-level rise on the order of 10’s of m/Ky (e.g. Fairbanks, 1989). These numbers are huge compared to other rates in the forward model originating from tectonic subsidence or tectono-eustatic. Glacio-eustatic events occurring with such rapid rates shall have a profound effect on accommodation space, and diagenesis.

Matthews (1988) compares the Haq et al. (1987) eustatic curves derived from sequence stratigraphy with those deduced from tectono-eustatic history and the deep sea oxygen isotope record (Figure 3). Matthews (1988) proposes the “long term” curve should be equated with ice-free world tectono-eustatic history and that the “short term” curve should be broken into a high-stand envelope and a low-stand envelope of the high-frequency glacio-eustatic signal. After identifying these differences in construct and terminology, the two approaches to sea-level history are in reasonably good agreement for the last 65 My. High sea-levels on the Haq et al. (1987) long term curve around 53, 34, and 15 million years before present (Ma) do not exist in Figure 4, nor are they discussed in any original sources cited by Haq et al. (1987) or Vail and Haq (1988). The major high-stand at 35-30 Ma in the sequence stratigraphy curves is weak to non-existent in the oxygen isotope record. Nevertheless, both constructs suggest a world intermittently ice-free prior to 45 Ma, a generally glacial world by 30 Ma, generally rising glacio-eustatic high-stand envelope 20-15 Ma, and generally falling glacio-eustatic high-stand envelope 15-6 Ma.

Calibrated in the Late Pleistocene, the Cenozoic deep sea oxygen isotope record becomes a proxy for sea-level history (Matthews and Poore, 1980; Prentice and Matthews, 1988, 1991). Further, it is encouraging that the sequence stratigraphy generalized curve is at least in partial agreement with the Cenozoic oxygen isotope record (Figure 3). Both of these records shall be used as calibration with regards to development of the concept of orbital forcing of glacio-eustasy.

Sea-level Estimates Based on Orbital Forcing Time Series

Orbital forcing of glacio-eustasy is an idea that has captured the attention of geologists for many years (Imbrie and Imbrie, 1979). Variation in seasonal distribution of solar insolation causes cyclic variation in global ice volume at frequencies derived from changes in the earth’s orbit about the sun. Predicted periodicities include precession (approximately 20 Ky); tilt (40 Ky); and eccentricity (100 Ky and 400 Ky) (Berger et al., 1992). The simplest model for glacio-eustatic sea-level variation as a function of orbital forcing comes from curve-matching among the deep sea oxygen isotope record, the late Pleistocene coral reef terrace record, and solar insolation time series at various latitudes (e.g. Broecker et al., 1968; Mesolella et al., 1969; Hays et al., 1976; Fairbanks and Matthews, 1978; Imbrie and Imbrie, 1980).

Analysis of Tertiary oxygen isotope time series suggest a different causal mechanism between orbital forcing and pre-Pleistocene Antarctic ice volume history (e.g. Prentice and Matthews, 1991). The formation of warm saline bottom water (WSBW) (Brass et al., 1982) in and/or around the Tethyan Seaway delivers sensible heat to high southern latitudes. Warm water upwelling near a cold Antarctic continent presents a great opportunity to grow large ice sheets. The greater the summer insolation at 30° North, the greater the sensible heat delivered to high southern latitude. Similarly, greater summer insolation at 70° South induces summer melting of ice accumulated on Antarctica during other seasons. This is the mechanism which most likely applies to the Mesozoic.

Figures 5 and 6 present the main features of our calculation of sea-level files from orbital forcing data. First, eccentricity, precession and tilt are calculated appropriate to the time interval under consideration (Matthews et al., 1997). Second, solar insolation (Figure 5, step 1) is calculated at critical latitudes and seasons. Third, a sea-level time series (Figure 5, step 2) is calculated with a dynamic model based on July insolation at 30° North (e.g. the WSBW production term) and January insolation at 70° South (the Antarctic summer ice melt term). Figure 6 provides examples of ice making/ice melting rates as a function of insolation at these two critical latitudes. Finally, results are compared to target sea-level proxy (deep sea oxygen isotope data or the sequence stratigraphy generalized curve) and model parameters are adjusted to achieve a best fit (Figure 5, step 3).

It is interesting to compare orbital forcing of glacio-eustasy with sequence stratigraphy sea-level generalities and nomenclature. The vast majority of “third-order cycles” are almost certainly the 2.0 My and 2.85 My cycles depicted in Figures 7 and 8. The 400 Ky eccentricity cycle is clearly displayed in the orbital forcing curve of Figure 7 and is recognized in practice in many sequence stratigraphy applications. It is clear in Figure 7 that each 400 Ky cycle contains higher frequency cycles. However, Figure 8 indicates no systematic relations among higher cycles within the 400 Ky cycles. The present authors shall continue to use “third-order” as synonymous with 2.0-2.85 My eccentricity nodes, and “higher-order” to denote 400 Ky and shorter period.

Choice of sea-level concept can have a profound effect upon static models derived from forward modeling. Figure 7 provides a comparison of classic sequence stratigraphy generalized curve and high-frequency glacio-eustasy calculated from orbital forcing. With regards to carbonate reservoir characterization and flow simulation, the major drawback here is that low-resolution concepts of sea-level fluctuation cannot reproduce realistic diagenetic overprint (e.g. Humphrey and Quinn, 1989). Indeed, the entire formalism of sequence stratigraphy is designed to recognize sedimentation units. Scarce mention is given to diagenetic overprint.

Thus, for forward modeling of carbonate stratigraphy, the sea-level input file of choice is orbital forcing of glacio-eustasy. Changes in boundary conditions may introduce long-term change to the general directions of the orbital forcing calculation. Aperiodic events may randomly disrupt the global ocean/atmosphere/cryosphere system. Nevertheless, through it all, there shall exist a high frequency signal which predicts both depositional stratigraphic sequence and diagenetic overprint.

Climate

Climate is the third of our top level independent variables (Figure 1). Temperature and moisture balance dictate the differences between skeletal grainstone, non-skeletal grainstone and evaporites. Meteoritic diagenesis shall require seasonal rains. Climate, like glacio-eustasy, varies with orbital forcing. Climate should also be expected to vary with the changing position of the continents back through time. These concepts are reviewed briefly below with special emphasis on the Arabian sub-continent and the Late Jurassic time interval.

Quaternary Climatic Variability

Large areas of the Arabian sub-continent appear to have been much wetter only a few thousand to tens of thousands of years ago (e.g. Roberts and Wright, 1993). Fresh water lakes existed in the Rub’ Al Khali where there are dry playas today. Lakes in excess of 20 m water depth with freshwater gastropods, ostracodes, and Chara are dated as older than 20 thousand years before present (Ka) by Carbon-14. Lakes 5 or 10 m deep occupy these localities intermittently around 9 to 6 Ka. Pollen data likewise suggests more rainfall and vegetative stabilization of sand dunes. Viewed from the perspective of orbital forcing of glacio-eustasy, these data suggest a significantly wetter climate for the Arabian sub-continent during the Last Glacial Maximum (LGM) and the transition from LGM to modern conditions.

Elsewhere, the Late Quaternary geology of Bonaire and Curacao, Netherlands Antilles, provides an example of modern evaporite sedimentation, side by side, with Late Pleistocene meteoric diagenesis, strongly suggestive of climatic variability on the time scale of orbital forcing. The islands are today quite dry. Both Curacao and Bonaire exhibit modern evaporite deposition up to the level of halite precipitation in restricted lagoons. Still, nearby subaerial Late Pleistocene constructive carbonate terraces (e.g. Herweijer and Focke, 1978) show abundant evidence of meteoritic diagenesis (Fairbanks and Mangion, 1974) including the development of karstic surfaces and sink holes (Murray, 1969). These data imply variation in moisture balance on orbital forcing time scales. Similar variation in the orbital forcing of climate might be expected in the Mesozoic of the Arabian Gulf.

Jurassic Paleoclimate of the Arabian Sub-Continent

Vastly different spatial configuration of the continents (e.g. Matthews, 1984, Figure 18.1) dictates that Jurassic climate may be different from today. Northern Arabia was near the equator and likely experienced seasonal rains from an intertropical convergence zone developed over a large ocean to the east. India was still part of Gondwanaland and situated at 30° South latitude. Thus, there was no Tibetan Plateau to generate today’s Indian monsoon.

Kutzbach and Gallimore (1989) presents a general atmospheric circulation simulation of Pangaean climates which is probably realistic for the climate of Arabia in Jurassic time. Precipitation exceeds evaporation for about half the year. Rainfall averages about 70 centimeters per year. Net evaporation is intense during June, July and August. On a seasonal basis, a net thickness of 2-4 m of water would evaporate off of shallow lagoons. On an annual basis, this number drops to a few tens of centimeters per year. While the results of this climate modeling experiment should be considered preliminary, they do illustrate the potential for using general circulation model output to quantify one-dimensional models of sedimentation and fresh-water diagenesis.

Sedimentation

While tectonics, eustasy, and climate are the major independent variables (Figure 1), it is computationally efficient to regard aspects of sedimentation and diagenesis as conditionally independent variables.

Regarding forward modeling of carbonate sedimentation, the problem conveniently breaks into two categories as follows: (1) what sediment types should be there in the first place (e.g. Schlager, 1981; Kendall and Schlager, 1981; Hallock, 1988; Hine et al., 1988); and (2) how fast can these sediments be produced under conditions which likely existed (e.g. Stockman et al., 1967; Neumann and Land, 1975). The types of carbonate sediment to be expected are often known from observational experience with the sequences under consideration. Recent analogs can serve to identify conditions likely to encourage the recurrence of various sediment types. Estimation of sedimentation rates is a more difficult problem.

What, Where, When, etc.?

Factors governing carbonate sedimentation are qualitatively well-known (Matthews, 1984). The warm, oligotrophic, western side of an ocean basin is the natural habitat of skeletal carbonates. At local scales, lime-secreting benthic organisms need substrate and room to grow, as well as protection from excesses of turbidity, nutrients, salinity, or temperature (Wantland and Pusey, 1975; Hallock and Schlager, 1986; Hallock, 1988). The carbonate platform interior may vary widely in environmental conditions and therewith sediment types. Given through-going circulation and a patchwork of shoal and channel bathymetry, skeletal carbonates may flourish (e.g. Southern Belize, Purdy et al., 1975). Given even seasonally anomalous salinity or temperature conditions, skeletal carbonates give way to oolite, pelloid and lime mud deposition (e.g. Great Bahama Bank, Purdy, 1963). Still further restriction leads to evaporite deposition (e.g. Butler, 1969; Kendall and Skipwith, 1969; Hardie and Shinn, 1986).

The above generalities constitute a snapshot of the generalized distribution of sedimentary facies at any one sea stand. When modeling stratigraphic sequences, the factors governing the accumulation of various sediment types are set in motion by rising or falling sea-level. It is convenient to think of a subset of sediment types as typifying a particular region within the paleogeography and then call upon water depth to determine which types occur in what order accompanying each sea-level change.

Sedimentation Rates

Regional studies in Belize (e.g. Wantland and Pusey, 1975) indicate that overall sedimentation rates are comparable for reef framework, transported reef flat detritus, and in situ subtidal skeletal carbonates. Carbon-14 dates from drill cores suggest that vertical accumulation rates on the order of 1 m/Ky are common whereas rates on the order of 3 m/Ky are seldom attained. Even these numbers are high compared to most estimates of sedimentation rate based on the study of ancient sequences. Is it possible that ancient sedimentation rates were consistently an order of magnitude lower for reasons not understood? Is it possible that modern and ancient sedimentation rates are comparable and 90% of geologic time is missing from the rock record?

Indeed, the question of missing time (e.g. subaerial exposure) and/or conditions of very low sedimentation rate (deposition of condensed section) complicate the task of estimating sedimentation rates for the forward modeling of carbonate sequence stratigraphy. Figure 9 summarizes a very simple hypothetical situation which illustrates the problem. Consider a very flat-lying carbonate paleogeography (shown in cross-section in the top panel) which is subtidal at all times on the left and intermittently subaerially exposed on the right. The lower panel depicts the occurrence of carbonate lithologies as they would occur within a single sea-level cycle at three paleogeographic locations. Note how the various sedimentary facies relate to sea-level history. To the left, the deep subtidal sediments occur at high-stand of the sea and the shallowest subtidal sediments occur at low-stand. To the right, deep subtidal sediments occur only during flooding of the previous subaerial exposure surface and are followed by shallow subtidal sediments, then by subaerial exposure. In the center, the deep subtidal sediments occur at high-stands, followed by shallow subtidal sediments as sea-level declines, then by subaerial exposure at extreme sea-level low-stand. These are all shallowing-upward cycles, but the relation between sediment types and sea-level history (and thereby chronostratigraphy and sedimentation rates) is different in each of these three paleogeographic locations.

Attacking the problem first from the landward side (the right portion of Figure 9), note that there are two ways in which time may not be represented in the resultant stratigraphic sequence. First, there is the actual absence of sediment deposition during times of subaerial exposure. Second, there is the matter of filling all accommodation space early in the transgressive phase. Commonly, if sedimentation shuts down, extensive cementation occurs at the sediment/water interface, resulting in the formation of a thin, dense hardground (e.g. Shinn, 1969). Such hardgrounds should be regarded as condensed sections, in that an unspecified amount of time is represented by cementation of existing sediment as opposed to continued deposition of new sediment. Both subaerial exposure surfaces and hardgrounds should be anticipated to have profound effects on permeability in relatively thin sedimentary layers.

At the basinward side of this problem (left side, Figure 9), hundred percent of time is represented by rock (i.e., no subaerial exposure is indicated). Nevertheless, there are at least two opportunities for condensed section to complicate the estimation of sedimentation rates. First, there is opportunity in shallow subtidal lithology for temporary conditions of non-deposition and the development of hardgrounds similar to those discussed above. Also, deposition of lime mudstone condensed section should be expected at the base of each transgressive phase of cycle development (e.g. Azer and Peebles, 1998). These units may be relatively thin and quite densely cemented, creating permeability barriers that might be missed on electric logs or in visual examination of core.

For the carbonate reservoir geologist, the problem is to decide which relationships depicted in Figure 9 operate within any particular stratigraphic section. Unfortunately, this distinction is commonly based upon electric log or visual examination of core. As noted below, it would be preferable to base this decision upon chemical and isotopic evidence of the diagenetic overprint produced by subaerial exposure, vadose and phreatic recrystallization, and mixed water dolomitization.

Diagenesis

Recent shallow-marine carbonate sediments consist primarily of the unstable minerals, aragonite and high magnesium calcite (e.g. Matthews, 1966; Husseini and Matthews, 1972), whereas carbonate rocks are primarily the stable minerals, low magnesium calcite and dolomite. Where and when mineralogical stabilization takes place can have a profound effect on the porosity and permeability of carbonate sediments and on the geometry of permeable and impermeable units within the carbonate reservoir. Here, early diagenesis in the presence of fresh-water is especially emphasized.

Rapid mineralogical stabilization, and related dissolution and cementation, occurs when shallow marine carbonate sediments are exposed to meteoric water (e.g. Matthews and Frohlich, 1987; Halley and Matthews, 1987; Budd and Vacher, 1991). As depicted in Figure 10, subaerial exposure of carbonate sediments within the geologic record is most commonly associated with glacio-eustatic rapid lowering of sea-level by a few meters to tens of meters. Subaerial exposure creates several important diagenetic environments. These are the subaerial exposure surface itself, and downward from it the vadose environment (pore space not saturated with water) and the phreatic environment (pore space 100% occupied by water). The phreatic can be further subdivided into the fresh-water, mixed-water, and salt-water phreatic environments.

The subaerial exposure surface itself is often characterized by development of a low magnesium calcite caliche profile (Harrison, 1974). The vadose environment is often the site of slow mineralogical stabilization. The fresh-water phreatic lense environment (when present) is the site of relatively rapid mineralogical stabilization and related dissolution and precipitation (e.g. Steinen and Matthews, 1973; Matthews, 1974; Allan and Matthews, 1977, 1982; Quinn and Matthews, 1990). The mixed-water phreatic environment is commonly the site of rapid dolomitization of carbonate sediments (e.g. Wigley and Plummer, 1976; Humphrey, 1988). This may occur where phreatic fresh-water mixes vertically or laterally with salt-water, or it may occur where vadose fresh-water mixes vertically with salt-water. The marine phreatic, away-from-the-sediment/water interface and away-from-mixing, is the site of extremely slow mineralogical stabilization (e.g. Major and Matthews, 1983).

From the standpoint of reservoir characterization, the vadose/phreatic interface can be the site of solution/precipitation reactions which result in pairs of strata exhibiting high permeability above, and low permeability below. This results from simple carbonate chemistry phenomena at or near the vadose/phreatic interface. So long as a volume of water resides in a mineralogically unstable vadose environment, calcite precipitation occurs relatively slowly. When this vadose water enters the phreatic environment, precipitation occurs and CO2 gas is released to the overlying vadose environment. The availability of still yet more CO2 gas, released from precipitation of carbonate mineral at the top of the phreatic environment, causes still more dissolution in the immediate overlying vadose environment. The result of this feedback is the development of micro-cavernous porosity on the scale of a few centimeters to meters immediately above the phreatic environment. A similar development of permeable/impermeable pairs can occur in the vadose environment where cementation of former subaerial exposure surfaces has created local aquicludes (Videtich and Matthews, 1980).

Reliable identification of diagenetic overprint and thereby geometry of permeable and impermeable units typically involves stable isotopes, minor element chemistry, and petrography. Forward modeling suggests the resulting carbonate sequences will exhibit a complex diagenetic history at the scale of a few centimeters to a few meters (e.g. Matthews and Frohlich, 1987). Several investigations have sampled carbonate stratigraphic sections at close sample interval (e.g. Allan and Matthews, 1982; Wagner and Matthews, 1982; Wagner, 1983; Humphrey et al., 1986; Halley and Matthews, 1987; Al-Sagri, 1989; Quinn and Matthews, 1990). These studies clearly indicate that diagnostic variation does indeed occur at the centimeter to meter scale. Further, there are no “magic numbers” with regards to diagenetic interpretation of stable isotope and minor element chemistry data. Rather, the pattern of stratigraphic variation holds the key. Indeed, diagnostic information may be contained even within the cement stratigraphy of individual rock samples (e.g. Benson and Matthews, 1971; Benson et al., 1972). Thus, there is no shortcut to data acquisition on individual cores. If a core is worth studying at all, it is worth studying at the scale of a few centimeters.

With regards to reservoir characterization, these observations pose a major question regarding sampling strategy. If “every core” is to be studied at the centimeter scale and with the full suite of analytical tools, the task becomes prohibitively expensive and time consuming. If every core is to be sampled at meter to five meter intervals, important stratigraphic variation shall surely go unsampled. The alternative is to adopt a sampling strategy based upon confidence in a process approach. Study a few key wells in great detail and use them as calibration for a forward model. Thus informed, make numerous forward model simulations and correlate throughout the field using more rapid, cost-efficient methods to map the geometry of permeable and impermeable units identified in the simulations.

Spatial and Temporal Considerations

Last, but not least, serious thought must be given to dimension and temporal precision of forward models. These important decisions control the number of calculations that must be made and the amount of information that must be stored during a model run.

Dimensions

Programs can be written to mimic geological processes in one, two, or three dimensions. One-dimensional models are particularly appropriate for carbonates because (1) sediment production commonly occurs in situ, (2) the sediment type and sedimentation rate can often be considered a function of very general paleogeography and water depth, and (3) the fresh-water/salt-water chemistry and hydrology controlling early diagenesis depends strongly on the single vertical dimension: elevation relative to sea-level. Thus, often the properties of a carbonate stratigraphic sequence can be inferred in detail from sea-level history without resorting to a two- or three-dimensional model (Matthews and Frohlich, 1987; Quinn and Matthews, 1990). Build to a two- or three-dimensional representation, one simulated stratigraphic section at a time. The goal is for the numerous one-dimensional models to provide detail and certainty as to what correlates from one spot to another and what pinches out somewhere in between. The goal is to postpone use of geostatistics to as far into the description of the reservoir as possible.

As a short cut to modeling carbonate reservoir properties, it may be computationally efficient to parameterize a one-dimensional carbonate model using results from a two- or three-dimensional model run. For example, a two-dimensional clastic model (e.g. Frohlich and Matthews, 1991) might be used to determine the distance from the shoreline to a shelf location where the carbonate reservoir facies occurs. This distance could then be parameterized as a proxy for local salinity and turbidity affects imposed by nearby clastic sedimentation, landward of the carbonate site. Similarly, three-dimensional global atmospheric circulation models (e.g. Barron et al., 1989; Kutzbach and Gallimore, 1989) might be used for determining moisture balance, seasonal temperature, and other important features of local climate.

Temporal Precision

An essential choice for the modeler is whether to distribute sediments at regular intervals of time (every one Ky), or only at sea-level still-stands. It is most efficient computationally to calculate sedimentation only at still-stands. Further, the highest high-stands and lowest low-stands often determine the locations of major stratigraphic boundaries. Alternatively, deposition of sediment at regular, equally-spaced time intervals is appropriate if a differential equation is used to distribute sediment. This is quite appropriate for clastic sedimentation and may fit some carbonate situations.

DISCUSSION

In this section, the generalities outlined above are brought to bear on the specific tasks of reservoir characterization and flow simulation in typical Arabian Gulf operating situations. The task is open-ended. A “history match” is a good beginning; but reservoir-performance forecasting is the goal. The science of reservoir characterization will likely continue to evolve. Construction of static models from dynamic forward models is a natural fit in this evolutionary process. Applications of these concepts to Arabian Gulf oil fields are discussed. The inverse problem is mentioned briefly.

The Natural Evolution of Reservoir Characterization

The science of reservoir characterization is approaching a transition which has occurred before in many fields of science (e.g., most of physics and chemistry). Early in the development of any scientific field, the emphasis is on observation of statics and the statistical summary of observed relationships. With experience, comes better understanding of the dynamics responsible for the observed statics. With understanding of dynamics, comes forward modeling.

So, where stands reservoir characterization and flow simulation now and where should the field be five or ten years from now? The typical procedure has been for geologists to construct a reservoir characterization static three-dimensional geocellular model (hereafter, static model). The rocks are regarded as objects to be understood qualitatively with regards to origin and described quantitatively with regards to porosity and permeability. Dolomite is commonly recognized as different from calcite. Calcites are commonly subdivided on the basis of grain size and grain type into such things as mudstones, wackestones, packstones, and grainstones (Dunham, 1962), organized under the rubric of sequence stratigraphy, and correlated throughout the field. Commonly, it is rock type in the context of sequence stratigraphy that is the initial characterization of each cell. Porosity and permeability are ascribed to the cells on the basis of statistical representation of rock types and observed spatial variability. Even the best rock types commonly exhibit wide variation in porosity and permeability. If no serious attention is given to processes responsible for this variability, it is common to distribute the good and the bad permeability more or less randomly within the cells assigned to each rock type.

There are some simple steps from static, stochastic modeling toward dynamic forward modeling. With processes and geological reasonableness in mind, construct alternative static models representing end members dominated by a particular process-related parameter (e.g., depositional facies, diagenetic overprint, fractures). Conduct sensitivity tests on several static models over the range of geological reasonableness. Several acceptable history matches will likely come out of this study. The questions now become whether future reservoir performance shall match only a subset of these static models and whether static models can be further improved by full implementation of dynamic forward models.

Construction of Static Models from Dynamic Forward Models

The forward modeling approach to the construction of static models offers room for improvement to the product, reduced costs, as well as opportunity to impose differing concepts in a systematic fashion. The main advantage is the explicit simulation of internal geometry of porosity and permeability imparted to the reservoir by depositional events and diagenetic overprint. For example, syn-depositional marine diagenesis commonly makes very thin, very tight hardgrounds. The effects of early subaerial exposure (Figure 10, location A) can be profound. The effects of early post-depositional freshwater diagenesis may improve permeability or may destroy permeability. These effects may be imparted to the sediments irregardless of original sediment type. The forward model can handle all of this from first principles, calibrated to the local situation via a few well-studied boreholes.

Geometry of depositional units and diagenetic overprint will likely need to be calibrated on a field-by-field basis. Nevertheless, some generic geometric relationships are self-evident. Marine hardgrounds and subaerial exposure surfaces follow paleotopography. Expect paleo-topography to rise and fall within the stratigraphic sequence. Expect some “para-sequence” units to correlate widely. Expect other units to have limited aerial extent. The tops of old water tables shall be flat at the time of their formation and shall correlate throughout the field irrespective of deposition-based sequence stratigraphy.

With field-wide correlation of sequence stratigraphy and diagenetic overprint geometry explicitly provided by forward model output, porosity and permeability can be parceled out within the static model with much greater explicitness. Tight rocks in grainstone facies almost certainly represent either submarine cemented hardgrounds, subaerial exposure surfaces, or tops of water tables; and so forward through the process of constructing the static model. In the process, geometric continuity of impermeable horizons and highly permeable intervals will replace the largely random distribution often assigned to cells by the descriptive/geostatistical approach.

Arabian Gulf Applications

Preliminary indications are that a process approach to reservoir characterization is applicable to Arabian Gulf carbonate reservoirs ranging in age from Permian through Eocene. In the cases outlined below, enthusiasm for this approach rests upon identification of profound early diagenetic alteration associated with high-frequency glacio-eustasy. The overriding theme is well-understood, but individual applications require local calibration on a field-by-field basis. The details of geometry and importance of various early diagenetic processes vary widely.

Application of a high density sampling program to a single well in Kirkuk field clearly identifies mixed-water dolomite as the major permeability-maker in the reservoir and top-of-phreatic-lense poikilitic calcite spare as a major permeability-spoiler (Al-Sagri, 1989; Halley and Matthews, 1987). In Kirkuk field, evidence of early meteoritic diagenesis is rampant. The upper portion of the carbonate section shows abundant evidence of vadose diagenesis and is substantially impermeable. Mixed water dolomitization obliterates original depositional fabric and creates permeability on the order of 200 millidarcies. Interbedded with mixed water dolomites are calcite grainstones which have escaped the dolomitization process. These are uniformly tightly cemented and impermeable. Interbedded with permeable dolomite horizons are dolomite horizons which subsequently have been reduced to low permeability by the deposition of poikilitic spare. The stratigraphic sequence of alternating permeable and impermeable lithologies extends over about 60 meters of section and is the major reservoir for Kirkuk field.

Interestingly, Herron et al. (1992) reports a very similar relationship of permeability to dolomitization for a Cretaceous study well in the Arabian Gulf. The results clearly indicate that dolomitization is responsible for the best permeability and that the alternation of calcite and dolomite is stratified in a manner similar to that observed in Kirkuk. Third-order and higher cycles appear to be well-recognized in the Shu’aiba Formation and are attributed to high-frequency changes in (glacio-?)eustatic sea-level position (e.g. Kendall et al., 1997). If third-order cycles are behaving nicely in the Arabian Gulf Cretaceous rocks, the higher-order cyclicity, which is at the heart of an orbital forcing approach to reservoir characterization, awaits elucidation regarding diagenetic processes operative in these sections.

Arab-D carbonates in Ghawar field also provide an interesting test of the ideas advanced in this paper. To begin, the three evaporites of the Arab Formation and the overlying Hith Evaporite Formation are excellent candidates for third-order cycles of sequence stratigraphy. As noted above, these cycles are well-represented in orbital forcing calculations of glacio-eustasy. If the third-order cycles impart such a strong depositional environment signal to this part of the section, might higher-order sequence stratigraphy cycles and diagenetic overprint geometry be present in Arab-D carbonates? Indeed, Bray (1997) cites subtle petrographic evidence of subaerial exposure throughout the Arab-D at Ghawar field. Further, preliminary calculations and comparison to published literature (e.g. Powers, 1962; Douglas et al., 1996; Hughes, 1996) are very encouraging. While evidence of widespread freshwater diagenetic overprint is lacking, massive dolomite units which commonly obliterate original depositional texture near the middle of the Arab-D carbonate unit (Douglas et al., 1996, Figure 9; Meyer et al., 1996, figure 11) are possible candidates for mixed water dolomitization during prolonged and complex third-order low-stand of sea-level. Differences in mineral solubility between mineralogically unstable vadose above and dolomite mixed-water phreatic below, create a positive feedback relationship across the top of the water table. There is precipitation of less soluble phases in the mixed-water phreatic environment below, and solution of more soluble phases in the vadose environment above. The results should be a zone of cavernous porosity/permeability on the scale of tens of centimeters to a few meters. Elsewhere in the section, permeable/impermeable pairs are expected where calcite cementation at old subaerial exposure surfaces creates local aquicludes within the vadose (e.g. Videtich and Matthews, 1980).

The Inverse Problem

This paper is laced with references to geological reasonableness and uniqueness. These problems are well-known to the solid-earth geophysicist under the general rubric of “inverse theory” (e.g. Menke, 1989). While a complete review of inverse theory is beyond the scope of this paper, brief comments below point in the general direction.

The complete task is (1) to devise a dynamic forward model based upon general principles applicable to the problem at hand, (2) to demonstrate that there exist model parameters such that a forward model solution is a satisfactory fit to observed data, and (3) to quantitatively evaluate the non-uniqueness of even the most “perfect” forward solutions. Items (1) and (2) have been dealt with in this paper. Item (3) begins with the simple idea of sensitivity testing among alternative static models, but rapidly becomes computationally challenging. Nevertheless, it must be recognized that this is a time of ongoing vast improvement to computational capabilities. Things that were major achievements five years ago will soon appear almost trivial.

These computational improvements have led to new methods for solving the inverse problem, even when the number of variables is too large for grid-search methods (i.e., sample multi-dimensional space for all possible forward solutions). The most promising methods involve “directed” random searches of the solution space, “genetic algorithm” methods and “very fast simulated annealing” methods (Sen and Stoffa, 1995). These begin with conventional Monte Carlo sampling of possible solutions and then apply direction in controlled steps to further delineate the range of solutions which are “nearby” to the solutions that most closely match the observations. In short, it is time to look beyond non-unique solutions to the forward problem and toward a planned approach to the solution of the inverse problem. Hopefully, this paper is a start in that direction.

CONCLUSIONS

The choice facing the oil production industry is whether to continue to pursue a data-intensive, labor-intensive, statistical approach to reservoir characterization or to replace it with greater reliance on understanding of processes responsible for permeable and impermeable units/horizons within the reservoir. With process, comes geometry. With geometry, comes less reliance on geostatistical representation of reservoir inhomogeneity; comes a better representation of reality.

The first step in deciding to try a process approach is to acknowledge that today’s final product has no claim to uniqueness. Color and 3-D do not impart uniqueness. A history match does not guarantee accurate simulation of future reservoir performance. Multiple, alternative static models and sensitivity testing are required to evaluate the uncertainty of reservoir-performance forecasting.

The great advantage of application of dynamic forward modeling technology to reservoir characterization is that the problem can be explicitly divided into numerous components which can be studied separately and quantitatively. Small changes in model input parameters can cause large changes in model output. Subsidence rate differences of one centimeter per thousand years can be the difference between a grainstone and a mudstone; a permeable sand and a tightly cemented hardground, etc. Numerous combinations of parameters will converge on target data and thereby form the basis for explicit alternative static models for sensitivity testing. These will be much more meaningful experiments than simply making arbitrary changes the geostatistical permeability algorithm.

The one-dimensional dynamic forward model is the modeling tool of choice. Calibrate the model on a few key wells. Contour the field with regards to key input parameters/target variables. Drill simulated wells at will to generate two- or three-dimensional arrays. What correlates from well-to-well and what does not, is known precisely because there is output of all relevant information for each depositional event in each simulation. It will also become apparent that a full calculation for every cell in a two- or three-dimensional array is computationally inefficient. Where correlation is simple, interpolation is sufficient. Where things are more complex, drill more simulated wells.

Finally, the increased efficiency of dynamic forward modeling should encourage reservoir-performance forecasting throughout the life of the field. Too often in the past, serious study of the reservoir was undertaken only when serious problems developed. With ever-increasing computing power and dynamic forward modeling, the goal should be to predict future problems and avoid them, rather than wait to observe problems and react to them.

ACKNOWLEDGEMENTS

The authors wish to thank D. Bice, G.C. Bond, H. Bosscher, T.A. Cross, A.W. Droxler, R.K. Goldhammer, P. Hallock-Muller, J.W. Harbaugh, M.T. Harris, C.G.A. Harrison, A. C. Hine, J.D. Humphrey, J.E. Joyce, C.G. St.C. Kendall, M.A. Kominz, D.T. Lawrence, M.A. Perlmutter, M.L. Prentice, T.M. Quinn, J.F. Read, D.J. Reynolds, W. Schlager, and D.L. Turcotte for correspondence and/or discussion regarding their modeling efforts and/or related research. We thank G.C. Bond, P. Hallock-Muller and M.A. Kominz for critical review of an early draft manuscript. We thank Moujahed Al-Husseini, Ibrahim Al-Jallal and four anonymous GeoArabia reviewers for helpful comments. While we have tried to give serious consideration to correspondents views and reviewers comments, we acknowledge that correspondents and reviewers do not necessarily share our views. Indeed, the senior author accepts full responsibility for any inaccuracies, biased opinions, ambiguities, and/or omissions. This work is based on research supported in part by the U.S. National Science Foundation and various oil companies.

ABOUT THE AUTHORS

Robley K. Matthews is Professor of Geological Sciences at Brown University, Rhode Island, USA, and is general partner of RKM & Associates. He has over 35 years experience in carbonate sedimentation and diagenesis and their application to petroleum exploration and reservoir characterization. Current interests center around the use of computer-based dynamic models in stratigraphic simulation.

Cliff Frohlich is Senior Research Scientist and Associate Director of the Institute for Geophysics, University of Texas at Austin and is a partner in RKM & Associates. He is best known for his research concerning deep-focus earthquakes. However, he also has an abiding interest in computer modeling of various mechanical and dynamical systems, ranging from bowling balls to stratigraphy of sedimentary basins.