Improved Reservoir Characterization in Low-Permeability Reservoirs with Geostatistical Models

SPE 28608

Improved Reservoir Characterization in Low-Permeability Reservoirs With Geostatistical Models

D.N. Meehan*, Union Pacific Resources Co, and S.K. Verma,* Stanford U.

* SPE Members

Copyright 1994, Society of Petroleum Engineers, Inc.

This paper was prepared for presentation at the SPE 69th Annual Technical Conference and Exhibition held in New Orleans, LA. USA, 25-28 September 1994.

This paper was selected for presentation by an SPE Program Committee following review of information contained in an abstract submitted by the author(s). Contents of the paper, as presented, have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material, as presented, does not necessarily reflect any position of the Society of Petroleum Engineers, its officers, or members. Papers presented at SPE meetings are subject to publication review by Editorial Committees of the Society of Petroleum Engineers. Electronic reproduction, distribution, or storage of any part of this paper for commercial purposes without the written consent of the Society of Petroleum Engineers is prohibited. Permission to reproduce in print is restricted to an abstract of not more than 300 words; illustrations may not be copied. The abstract must contain conspicuous acknowledgment of where and by whom the paper was presented. Write Librarian, SPE, P.O. Box 833836, Richardson, TX 75083-3836, U.S.A., fax 01-972-952-9435.

Introduction

Abstract

Infill drilling has been commercially successful in many low permeability, heterogeneous gas reservoirs. Reservoir discontinuities have often been suspected as a factor in poor gas recoveries on wide spacing. Large vertical and lateral variations in permeability make it difficult to account for partial drainage at infill locations. How many wells must be drilled to recover the gas? What are the effects of heterogeneities on optimal well spacing and fracture length?

In this paper, a case history illustrates the power of incorporating high resolution, fine grid geostatistical models in simulating reservoir behavior. Previous reservoir simulation studies provided acceptable matches of flow rates and pressures by fairly arbitrary reductions in the log derived net pay for the entire reservoir or away from the well. However, these models failed to match extended pressure buildups. The buildups indicate significantly higher gas-in-place in the reservoir than is indicated by simulation matches based on simpler reservoir descriptions. The geostatistical model presented here resulted in excellent multiple well history matches and matched the long term buildup.

The techniques for generating the reservoir description are summarized along with the reservoir simulation results. Predictions of infill drilling success with this model are better than for prior models. Predictions of incremental gas recoveries from infill drilling from this model are consistent with observed results. Reservoir heterogeneities (specifically the lateral continuity of permeability) appear to be the most important factors in this reservoir controlling inadequate drainage of the uppermost intervals. These lateral heterogeneities appear to be diagenetic permeability alterations resulting in partial compartmentalization of the many individual sands.

Optimal well spacing in very low permeability reservoirs has been addressed by numerous authors. Wells with permeabilities in the Cotton Valley range (# 0.01 md) generally indicate extremely long "optimal" fracture lengths (often in excess of 1000 ft fracture half-lengths). It is doubtful that such fractures can be created without vastly larger jobs than predicted by conventional hydraulic fracture models. Economic approaches used in conventional fracture optimization models may be inappropriate in thick intervals with few stress barriers. Inadequate barriers to fracture height growth and reservoir heterogeneities indicate the need for closer spacing and moderate fracture lengths.

Continued infill drilling accomplishes two things, viz. increased access to poorly drained portions of the reservoir with better stimulations and acceleration of recoveries from the most continuous portions of the reservoir. Current well costs can justify incremental recoveries at the current spacing levels; however, significant gas will remain unrecovered. The importance of lowering well costs is described.

Geology

The Carthage (Cotton Valley) Field is located in Panola County in East Texas. The Cotton Valley sandstones of Carthage field consist of a series of marine and lagoonal deposits overlying the gentle regional structure associated with the Sabine uplift. At its crest the Cotton Valley section is 1200 ft thick, expanding to 1500 ft downdip. The Cotton Valley interval includes very fine-grained sandstones, siltstones, shales and limestones. Sediments were deposited by longshore currents that deposited continuous clean sands in a shallow marine environment. Shale laminations are extensive, resulting in small sand members ranging in thickness from a few inches to 10 to 15 feet. Bounding shale laminae are lenticular and discontinuous. Diagenesis in the form of calcite cementation and quartz overgrowth, combined with overburden pressure have dramatically reduced porosity and permeability. Sand porosities range from 2% to 12% with microdarcy-level permeabilities. Massive hydraulic fractures stimulations are required for commercial completions.

Previous Studies

Modern hydraulic fracturing techniques and improved natural gas prices resulted in rapid development of the Carthage (Cotton Valley) gas field in the 1976-1979 time frame. Attempts to model well performance followed quickly, with well test and simulation studies indicating hydraulic fracture lengths much shorter than predicted by conventional 2-D fracture models. One reservoir simulation study was completed in 1992 by an in-house team on the Carthage Gas Unit 21 (CGU-21) to evaluate 80-acre drilling potential. Cartesian grids with one of the directions oriented in the expected fracture direction were used with uniform reservoir properties. The reservoir was divided into three non-communicating layers and one communicating layer. A history-match with well head pressure controls was performed. The only way that a good match could be obtained was by compartmentalizing the upper layers. A similar approach was taken by Meehan and Pennington1 and by Schell2. Individual flowing pressure declines were matched; however, the model pressure could not increase to the observed value in well 21-2 when it was shut in for a eighteen month pressure buildup. Simple single layer models were also made but required large decreases in net pay and decreasing reservoir permeability with time (or increasing skin effect) to match production.

This Study

The study was divided into four stages,

1. Data analysis,.

2. Reservoir characterization based on geostatistical methods,

3. Creating a reservoir model for numerical simulation, and

4. Matching reservoir and well performance and making reservoir performance predictions.

Data Analysis

Two types of data were used in the study: data to arrive at the geological model of the reservoir and production/pressure data for each well. The first stage of the project consisted of gathering all log, core, production and pressure data. Logs were recalibrated and interpreted on a consistent basis, matching core data for porosities and shale content. Flowmeter logs were available at several times for most wells; individual layer flow rates were used as history match parameters. Flowing tubing pressures, well tests, pressure buildups, pressure gradient checks and production data were also integrated. Well flow rates and flowing tubing pressures were used to calculate bottom hole flowing pressures. Incorporating flowing gradient data improved the pressure drop calculations.

Geological Data

Exhaustive petrophysical studies of all wells incorporating the full range of available open-hole logs and core analyses were conducted. Foot-by-foot estimates of porosity (N), shale volume (Vsh), and water saturation (Sw) were made and integrated with formation tops and bottom.

Analytical plots (histogram, probability, scatter, etc.) for N , Vsh and Sw were made for each group and sub-group of sands. This analysis is useful to understand frequency distributions, detect correlations between properties, identify outliers, provide regional statistics, etc. There is only a modest correlation between porosity and water saturation for most groups of sands. These plots were also useful in preparing geostatistical simulations to be undertaken and in understanding the numerous realizations.

Production and Pressure Data

The section taken for study was around CGU Unit 21. The area covered by the study is 9000 ft in the X-direction and 7000 ft in the Y-direction (1446 acres). There are ten wells in this area which have produced from the Lower Cotton Valley Sands. The surrounding 24 wells were not included in the reservoir simulation match but were analyzed for variogram development and geostatistical modeling. The first well (21-2) in the simulation area commenced production in January, 1979. This well produced intermittently due to gas demand. Early bottom hole pressure buildup measurements failed to stabilize. However, wellhead flowing pressures were available for entire well history along with numerous measured bottom-hole pressures. Measured and modeled bottom-hole pressures were in excellent agreement, providing confidence in using flowing tubing pressure in the simulation runs.

Geostatistical Simulations

The geological model was generated using a geostatistical approach provided by a group at Stanford University based on the "Amoco Data Set" 3. This approach is summarized here without extensive discussion of the geostatistical principles involved. The first step was to get facies distributions, followed by determination of N, Vsh and Sw. This information was used to provide initial estimates of permeability. Formation tops for each interval were determined by kriging.

Vsh as a function of areal location and depth was the first attribute addressed. Following our statistical study, we examined the spatial continuity of each reservoir property as measured by the variogram. Variograms are a first-order measure of an attribute's spatial variability.

Computing the variogram

Spatial variability is commonly measured by the semi-variogram, defined as the average squared difference between two attribute values approximately separated by vector h :

where N(h) is the number of pairs, xi is value at start or tail of pair i and yi is the corresponding end or head value. h can be specified with directional and distance tolerances. A semivariogram is normally used for the same variable, e. g. two N values separated by h.

Another useful measure of spatial variability is the indicator semivariogram. This variogram is computed on an internally constructed variable and requires the specification of a continuous variable and cutoff to create an indicator transform. For a specific cutoff and datum value the indicator transform is defined as:

Horizontal and vertical indicator semivariograms of Vsh for each group of sands was computed. Cutoffs were based on the cumulative probability distribution of the variable. Widely available GSLIB programs2 were used to compute the variograms as well as perform all the geostatistical modeling used in this study.

Modeling the Variogram

Standard variogram models easily fit the data; example computed and modeled horizontal indicator variograms show the vertical variogram of Vsh for one group of sands (Figure 1). All the wells in the simulation study area were used to develop the variogram models along with the offset wells within 3000 ft of the simulation area. For most horizontal variograms a spherical model was sufficient to model horizontal variability while a combination of exponential and spherical was used to model vertical variability. Not all cutoff levels show good horizontal correlation; vertical variograms are better correlated because of the presence of short-scale data. Data in the x-y plane are sparsely located with the minimum distance between two wells being about 900 feet. At some cutoff levels a model with range greater that 900 was observed. For levels where the correlation range from the available data was not observed, a value less than 900 ft was used.

Vsh Estimation

An estimate of a property (V)at any particular point can be made by a linear combination of values of the property at a set of given data points(Vi).

The challenge is in finding best possible weighting factors (wi) to be used with available data. One method is to assume a stationary random function as our model and specify its variogram. Taking this model as a true representation the values of which minimize the error variance are used to find V. Thus the variable to be minimized is:

where ri is the error of the i-th estimate and mr is the average error. If all the available data points are used at once then one does ordinary kriging. Ordinary kriging provides the best linear unbiased estimate and gives a very smooth picture and is in fact a contouring technique. Kriging provides a single numerical model which may be considered best in a local accuracy sense.

Stochastic simulation, on the other hand, is the process of drawing alternative, equally probable, joint realizations of a variable from a random function model. The realizations represent a number of possible images of the spatial distribution of the attribute values over the field. Each realization, also called a stochastic image, reflects properties that have been imposed on a random function model. Typically, the realizations honor input attribute values at data locations and are thus said to be conditional. Such conditional simulations correct the smoothing effect shown on maps produced by the kriging algorithm. In the sequential simulation approach all original data in a given neighborhood (of the point where the property is to be estimated) as well as all simulated values available up to that point in the simulation are used to obtain the estimate. The size of the conditioning data set increases as the number of values at simulated points increases. As far as the implementation of such algorithms are concerned, a limit is set on number of original and simulated data that can be used to obtain the estimate. A random sequence is followed in selecting the nodes where the property is to be simulated. When indicator semivariograms are used for the random function model then the process is called sequential indicator simulation5.

Numerous realizations were obtained by changing the random path followed in the simulations. Figure 2 gives three such simulations of Vsh for one group of sands. Similar simulations were performed for each of the other five major groups of sands. Color output is essential to properly visualize these results.

Indicator simulations were performed on a grid of 45 by 35 blocks in the x-y plane. These grids were 200 ft in length in each direction. In the vertical direction, the stratigraphic thickness of each group of sands was used. Grids in the vertical direction were five feet thick.

Selection of Realization

Any of the realizations shown in Figure 2 could be selected. Selection of the best possible realization may be very important. In the current study, available data were distributed throughout the area of interest without excessive local concentration. Hence it was decided to select that realization which had a similar cumulative probability density function as the data. Quantile-Quantile (q-q) plots (Figure 3) and histogram plots for each of the realizations was made. The realization which gave a best fit around the 45 degree straight line on the q-q plot was selected. A similar approach was used for each group of sands.

Modeling Porosity

It was observed that for all the groups of sands there was a good correlation between Vsh and and a good crossvariogram existed (at least in the vertical direction). Porosity realizations were initially modeled independent of the Vsh realizations. This approach did not result in an acceptable correlation between the realizations for Vsh and N. Following the approach outlined in Ref. 5, it was necessary to make the N realizations dependent on the Vsh realizations. Markov-Bayes6 simulations were used for N to account for the relationship between Vsh and N by using Vsh data as soft indicator data. This approximates indicator cokriging, where the soft indicator covariances and cross-covariances are calibrated from the hard indicator covariance models.

Modeling Sw, Permeability and Formation Tops

Sw values at each grid location were originally generated using Markov-Bayes simulations with the Vsh and values as soft data points; however, linear correlations of Sw with and Vsh were found to reduce computational time and generate very similar results. Vsh, N, and Sw values were used to determine an initial permeability value for each grid point using a relationship of the form:

where  is a constant obtained by history matching. Permeability values were modified during the history-match; a typical realization for permeability is given in Figure 4. Formation tops for individual layers were obtained using ordinary kriging.

Numerical Reservoir Modeling

Fine scale realizations of Vsh ,Sw , N, and permeability were computed at grid nodes whose dimensions were 200 by 200 by 5 feet (385,875 grid points). Flow simulation grid point locations had to be reduced to solve the problem on a workstation in a reasonable amount of time.

Upscaling

In the horizontal plane a decision was made to stay with the 200 ft by 200 ft block dimensions (the scale at which geostatistical simulation was done) to keep enough blocks between infill wells.

Simulator performance was determined to be acceptable with up to 40,000 grid blocks (a few hours per run), dictating the level of vertical upscaling. Layers were grouped to lump high Vsh content (shaly) intervals reducing the model to 24 layers with 37,800 grid blocks.

Simple upscaling techniques were used for computing effective permeability of the coarse blocks because the reduction factor was only about 0.10 and adjacent fine layers of similar Vsh properties were grouped together. Vertical permeability was computed by harmonic averaging with horizontal permeability computed using arithmetic averaging. Effective porosity of the coarse blocks was also computed by arithmetic averaging.