Effects of Shales in Fluvial-Deltaic Deposits: Ground ...

3 downloads 0 Views 1MB Size Report
Djuro Novakovic2, Christopher D. White2, Rucsandra M. Corbeanu3, William S. ...... Combination of arithmetic and harmonic averages (Li, Beckner and Kumar, ...
Mathematical Geology

Ms. No. ___-___

Effects of Shales in Fluvial-Deltaic Deposits: Ground-Penetrating Radar, Outcrop Observations, Geostatistics and Three-Dimensional Flow Modeling for the Ferron Sandstone, Utah1 Djuro Novakovic2, Christopher D. White2, Rucsandra M. Corbeanu3, William S. Hammon III3, Janok P. Bhattacharya3, and George A. McMechan3

ABSTRACT Ground-penetrating radar (GPR) surveys, outcrop measurements, and cores provide a high resolution threedimensional geologic model to investigate the effects of shales in marine-influenced lower delta-plain distributary channel deposits within the Cretaceous-age Ferron sandstone at Corbula Gulch in central Utah, USA. Shale statistics are computed from outcrop observations. Histograms of shale length and variograms of shale occurrence on accretionary surfaces indicate only slight anisotropy in shale dimensions. The estimated mean shale lengths in paleoflow-parallel and paleoflow-perpendicular directions are 4.9 and 6.0 m, respectively. The anisotropy in mean shale length (L/W = 1.2) is similar to the ratio of the integral ranges of Gaussian correlograms for the shale indicator in the corresponding directions (L/W = 1.3). Truncated Gaussian simulation generates maps of shales that are placed on variably dipping stratigraphic surfaces interpreted from high-resolution 3D GPR surveys, outcrop interpretations and boreholes. Sandstone permeability is estimated from radar responses calibrated to permeability measurements from core samples. Experimental design and flow simulations examine the effects of variogram range, shale fraction, and trends in shale fraction on predicted upscaled permeability, breakthrough time, and sweep efficiency. Approximately 1500 flow simulations examine three different geologic models, flow in all three coordinate directions, 16 geostatistical parameter combinations, and 10 realizations for each combination of parameters. Analysis of the flow simulations demonstrates that shales decrease sweep, recovery, and permeability, especially in the vertical direction. The effect on upscaled horizontal flow properties is less significant. Flow predictions for ideal tracer displacements at Corbula Gulch are sensitive to shale fraction, but are not sensitive to twofold variations in variogram range or to vertical trends in shale fraction. Keywords: distributary channel, reservoir simulation, conditional simulation, reservoir analog

1

Submitted for review August 27, 2001. The Craft and Hawkins Department of Petroleum Engineering, Louisiana State University, Baton Rouge, LA 70803 3 Department of Geosciences, The University of Texas at Dallas, P. O. Box 830688, Richardson, TX 75083-0688 2

2

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

INTRODUCTION Problem Statement This study examines the effects of shale drapes in laterally accreting point bar deposits. Many investigators have recognized the importance of low-permeability shales in controlling flow behavior. Shales affect vertical permeability (Begg and King, 1987), sweep efficiency and breakthrough time (Jackson and Muggeridge, 2001), and upscaled multiphase flow properties (Narayanan, 1999). However, it is difficult to include realistic distributions of shales in reservoir models because of the wide spacing of wells and limited vertical resolution of seismic surveys. One solution is to characterize the shale distribution statistically based on data from outcrop exposures of analogous deposits. Desbarats (1987) computed indicator variograms based on shale maps from an outcrop exposure of the Assakao sandstone for use in analogous fluvial reservoirs; these calculations assumed uniform shale dip. Geehan and Underwood (1993), Visser and Chessa (2000), and White and Willis (2000) present methods to compute unbiased estimates of shale length histograms. Because they are object-based, meaningful length distributions can be calculated even if the depositional dip varies. However, none of these studies addresses shale placement in a model with complex geometry such as point bar deposits, where the shales are deposited along surfaces that have variable dip. In addition, many of these studies were confined to two spatial dimensions (horizontal and vertical). In clastic reservoirs, shales are often not horizontal. Depositional dip may vary areally, and shale dimensions may be anisotropic. Methods are needed to create geologic models with such shale geometries and to predict their effects on flow behavior.

Approach A stratigraphic framework specifies the surfaces on which shales might occur. Shales are modeled geostatistically in two dimensions on the surfaces. Statistics including length distributions, shale indicator variograms, and shale proportions are computed in a stratigraphic coordinate system. Trends in shale proportion are also quantified along and between surfaces. In this study, the stratal surfaces are imaged by a high-resolution ground-penetrating radar survey (Corbeanu, 2001). These are the surfaces on which shale statistics are inferred. Corbeanu (2001) presents variograms of radar amplitude on stratigraphic surfaces. The 3-D GPR survey used nominal 100-MHz data, which has a depth of penetration of approximately 10-m and vertical resolution of a few decimeters (Corbeanu, 2001). The

Effects of Shales in Fluvial-Deltaic Sediments

3

survey grid is positioned to subcentimeter accuracy using global positioning system technology (Xu, 2000). Because lithology is correlated with radar responses (Hammon and others, 2001), variograms based on radar responses are related to the variograms for variables such as the shale indicator. However, no physics-based analysis of radar responses has linked the spatial correlation of radar amplitude to the spatial correlation of shale occurrence because the requisite data are unavailable. The errors of using amplitude as a proxy for shale occurrence have not been assessed. The errors due to imperfect correlation of amplitude with shale occurrence will add to the errors in variogram inference and stochastic simulation (Deutsch, 1997). Furthermore, the thin shales on the accretionary surfaces (less than the decimeter resolution of the 100-Mhz GPR data) are not resolved by ground-penetrating radar data. Therefore, variograms for the occurrence of thin, bed-draping shales are computed from outcrop exposures in this study. Although the geostatistical model for shale drape occurrence is based on observations from outcrop exposures, there is uncertainty in factors such as variogram range. The effects of these factors are addressed via an experimental design in which four factors (strike or paleoflow-perpendicular range, dip or paleoflow-parallel range, shale drape fraction, shale fraction trend) are varied over feasible ranges. The effects and interactions of the four factors are quantified with a four-factor two-level full factorial design (Box et al., 1978). Because the indicator models for shale are stochastic, multiple realizations are considered. Ten realizations are created for each factor combination in this study. Geologic models based on GPR surveys of two different volumes are analyzed, and flow is simulated for all models in three orthogonal directions. The flow simulations use a stratigraphic cornerpoint grid (King and Mansfield, 2000) and a commercial reservoir simulator (Schlumberger Technology Co., 1997) to predict the behavior of an ideal tracer displacement. Several flow responses are considered for each flow simulation: the upscaled permeability, based on a steady-state pressure solution; the time at which the produced tracer concentration exceeds 1 percent (often called breakthrough); and the fraction of the model contacted by injected fluid after 1 pore volume of injection (sweep efficiency). These flow responses are examined with linear sensitivity analysis of the two-level factorial (Myers and Montgomery, 1995). This approach ranks the effects and characterizes interactions between factors statistically (e.g., Kjønsvik and others, 1992). This analysis can be used to compute effective properties (Narayanan, 1999) and to estimate the value of more accurate characterization of shales (White and others, 2000).

4

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

Geologic Setting The statistical and flow models in this article are based on outcrop, core, borehole, and ground penetrating radar studies of the Ferron sandstone, which is a lithostratigraphically defined member of the Mancos Shale Formation. The Ferron sandstone was deposited in the Western Interior Seaway during Late Cretaceous (Turonian) time (Armstrong, 1968). The Ferron sandstone contains the preserved deposits of rivers, deltas, lagoons, shorefaces and associated deposits including coals. It is divided into seven sandstone-rich tongues, each overlain by a coal (Fig. 1; Ryer; 1983; Gardner, 1995). The Ferron sandstone has been used as an analog for fluvial-deltaic reservoirs, especially for the river-dominated systems of the U.S. Gulf of Mexico including the Frio sandstone (Barton, 1994). The study site is located south of Interstate 70 within the outcrop belt of the Ferron sandstone along the western margin of the San Rafael Swell in central Utah, USA (Fig. 2). The study location, Corbula Gulch, is interpreted to be marine-influenced point-bar deposits of distributary channels within the lower delta plain (Corbeanu, 2001). The deposits at Corbula Gulch include trough-cross-bedded and ripple-cross-laminated sandstones, sandstones, mudstones, shale clast lags, and thin mud drapes. The Corbula Gulch study site is within stratigraphic cycle 7 of Gardner (1995; Fig. 1). At the time of deposition, the delta shoreline was 15 to 20-km northeast of the study locality (Fig. 2). The heterogeneities most likely to affect flow in this system are mudstone bodies and thin, discontinuous mudstone drapes and mud-chip lags associated with accretionary surfaces of the point bar deposits. The combination of outcrop exposures and the high-resolution stratigraphic framework obtained from interpretation of groundpenetrating radar studies allows systematic investigation of the effects of these surfaces in three dimensions.

GEOSPATIAL MODEL Stratigraphic Framework Ground penetrating radar can image geologic features to depths of about 10-m with resolution from a few decimeters to approximately one meter, depending on the characteristic frequency (Corbeanu, 2001). Several surveys were performed at the Corbula Gulch locality; two fine-grid surveys, labeled A and B, are the subjects of this study (Fig. 3). Their dimensions are 27 by 31m and 51.5 by 28 m, respectively (Table 1). The acquisition and interpretation of these surveys have been discussed previously (Corbeanu, 2001; Hammon and others, 2001).

Effects of Shales in Fluvial-Deltaic Sediments

5

The interpretation of the GPR and outcrop data sets includes identification and description of a set of relatively continuous surfaces. The shape of these surfaces and correlation with surfaces exposed in outcrops supports their interpretation as accretionary surfaces within the point bar deposit (Corbeanu, 2001). Six surfaces were interpreted in the volume beneath grid A, and ten in the volume beneath grid B. These surfaces were estimated on GPR data volumes with horizontal spacing of 0.50 m. For flow modeling these surfaces are resampled onto grids with a spacing of approximately 1 m horizontally. Resampling reduces the gridded flow models to sizes that are feasible for flow simulation of multiple realizations and geologic models. The stratal frameworks for both grids have the inclined surfaces characteristic of point bar deposits (Fig. 4). The longer (31-m) side of grid A is approximately perpendicular to paleoflow (compare Figs. 3, 4A), and this is reflected in the offlapping geometry. In comparison, grid B (Fig. 4B) is oblique to paleoflow, and the dipping surfaces therefore appear more complex. The different orientations of the two surveys permit examination of flow responses in paleoflow-parallel, paleoflowperpendicular, and two oblique directions.

Estimation of Rock Properties To construct a flow model, rock properties such as porosity and permeability must be assigned. Flow properties are estimated by correlating them to radar responses (Hammon and others, 2001; Corbeanu, 2001). Permeability is correlated separately for each depositional unit. Within each unit, the instantaneous amplitude and phase of the radar responses are clustered by facies based on cores and borehole GPR data. The clusters are identified from the borehole data and are used to assign facies and permeability within the GPR data volume. The logarithm of fluid permeability is estimated as a linear function of instantaneous phase and frequency (Corbeanu, 2001). Permeability images are shown in Fig. 5. Porosity, which does not vary as widely, is estimated by correlation with the estimated fluid permeability based on core plug data from on-site boreholes and outcrops. The variability of estimated sandstone permeability is not large, covering only about a factor of two. The small range is partly due to the linear estimation method and resampling. The relatively high permeabilities toward the bottom of both grids will strongly affect flow behavior. Some of the permeability models used in this study include mudstone bodies, whereas other models do not distinguish regions of mudstone from surrounding sandstone (Corbeanu, 2001). Where distinguished, the mudstone bodies were identified from radar responses (previous paragraph). The mudstone bodies at Corbula Gulch average

6

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

greater than 1 dm in thickness, and therefore can be detected with GPR. The lateral extent of these bodies is usually less than 3 to 5 m, and mudstone bodies may be either within sandstone bodies or on accretion surfaces. Mudstone bodies may be true shales or mud-chip lags; these facies are difficult to distinguish with GPR data. In contrast, shale drapes on accretionary surfaces are often too thin to be interpreted from GPR data (less than 1 dm), but these drapes are sometimes laterally extensive. Throughout this article, “mudstone body” refers to a shale that is resolved from GPR data, and “shale drape” refers to a shale that that is below the resolution of GPR data and that lies on an accretionary surface. Because the shale drapes cannot always be detected using GPR data, they are modeled stochastically (Geostatistical Models of Shale Drapes, below). Some shales on accretionary surfaces are thick enough to resolve with GPR data; these thicker shale drapes are included in models as deterministic mudstone bodies. Placing a stochastic shale drape on a surface that is already covered by a mudstone body is redundant but will not lead to errors. It is consistent with the outcrop statistics and GPR data, and flow is blocked whether a shale drape, a mudstone body, or both a drape and body lie at a particular location on an accretionary surface. The fluid permeabilities of the mudstone bodies and the shale drapes are assumed to be zero. Flow models for grid A are constructed with both permeability models (including and not including mudstone bodies). If mudstone bodies are included, they fill 7.6 percent of grid A. The mudstone bodies appear as blank regions in Figure 5C. For grid B, only the model neglecting mudstone bodies is considered. In models that do not include mudstone bodies, the sandstone is more continuous because no zero-permeability mudstone bodies are incorporated in the grid (cf. Fig. 5A to 5C). The stratigraphic frameworks of the drapes-only and mudstone models are identical (Fig. 4); only the rock properties change. Including shale drapes does not alter the stratigraphy because they are modeled as infinitely thin barriers on accretionary surfaces. Computationally, the shale drapes are included in the flow models as z-direction transmissibility multipliers that prevent flow in the z-direction but do not alter the x- or y-direction properties of adjacent gridblocks.

Grid Construction for the Flow Model The interpreted surfaces from the GPR data volume are edited with stratigraphic modeling software. A cornerpoint grid specification (Ponting, 1992) is used to ensure that the stratigraphy is represented accurately, which cannot be done with conventional block-centered grids (King and Mansfield, 1999). Use of conventional grids for inclined geologic features may be inaccurate (Willis and White, 2000). Six and ten interpreted surfaces define the

Effects of Shales in Fluvial-Deltaic Sediments

7

stratigraphic framework for grids A and B, respectively (Fig. 4). The framework surfaces are edited to extend over the entire grid. Gaps are filled by linearly extrapolating beyond the input surfaces. The intersections of surfaces define pinchout locations, and surfaces are smoothed to avoid problems computing transmissibilities for very irregularly shaped blocks (King and Mansfield, 2000). Additional surfaces are defined to interpolate the stratal geometry in a geologically reasonable way, to capture permeability heterogeneity more accurately, and to improve flow model accuracy. The infilling surfaces are constructed parallel to the top bounding surface to mimic the accretionary geometry of these deposits. The geocellular grids have from 218,072 to 615,752 active gridblocks. Because these gridblock counts are too large for multiple reservoir simulation runs to be feasible, the grids are resampled at intervals of 1 m horizontally and approximately twice as coarse as the geocellular model vertically. Statistics of the resampled grid are given in Table 1 (48,960 to 93,630 active gridblocks). The resampling is kept as uniform as possible to avoid introducing differences into the models that reflect discretization rather than geology or flow behavior. The coarsening factor is not exactly eight (23) because the vertical grid is an irregular cornerpoint grid. The resampled grid sizes are feasible for numerical simulation; these are still large grids (Nt > 104) but are manageable (Nt < 105). Simulator time (or 3

work) requirements increase approximately proportionally to

N t 2 , where Nt is the number of active blocks.

Reducing the resolution by a factor of eight changes model run time by a factor of 23, enough to make a formerly day-long run complete in just over an hour. These savings are useful if many flow models are being simulated. The regularly spaced fluid permeability arrays derived from radar properties are resampled to the centroids of the cornerpoint gridblocks (Fig. 5). The grids capture the inclined features observed in the outcrop and in the GPR surveys. This level of resampling maintains much higher resolution than that typically used in reservoir simulation studies. The mean block volume for these models is about 0.25 m3 (Table 1). One relatively small gridblock in a field scale simulation model would have a volume of approximately 1000 m3 , about one-tenth the volume of the entire model for grid A (48,960 gridblocks). The large disparity between scales for geologic models and field scale flow models is one of the principal motivations for upscaling.

8

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

GEOSTATISTICAL MODELS OF SHALE DRAPES Indicator (Desbarats, 1987) and object (Haldorsen and Lake, 1984; Begg and King, 1987) models have been used to describe and to model shales. Object models use histograms of object dimensions, orientations, and locations; these statistics may or may not be correlated. The indicator approach is based on two-point covariances. The statistics discussed below are based on measurements from two cliff faces (Fig. 3). The N-S cliff face is nearly parallel to paleoflow. The E-W face is nearer to perpendicular to paleoflow. The E-W data are analyzed as if they were perpendicular to paleoflow, because no data are available exactly perpendicular to paleoflow. The paleoflow-parallel and paleoflow-perpendicular directions are assumed to be the principal correlation directions for this system.

Histograms Because the shale drapes are too thin to be resolved in the GPR data, drape lengths are tabulated from highresolution photomosaics of the two cliff faces (Fig. 3). Along each major stratigraphic surface, the presence or absence of shale is recorded at 1 m intervals. These observations are used to create histograms of shale length (Figure 6A) along the two cliffs. Fifty thin shales (or fragments of shales) occur along the accretionary bedding planes in the outcrop exposures. The mean observed shale lengths are 4.3 and 5.6 m for the N-S and E-W cliffs. There is no pronounced anisotropy in the histograms, and the termination frequencies t for the N-S and E-W faces are similar (0.41 and 0.33 m-1, respectively). The mean length, m L , can be estimated from the termination frequency, t, using the formula

mL =

2 . This estimate of the mean corrects for bias in length observations on an t

irregularly shaped outcrop exposure (White and Willis, 2000). This yields unbiased estimated means of 4.9 and 6.0 m on the N-S and E-W cliffs, respectively. The slightly greater mean length estimated from the termination frequency models is caused by a correction for partial exposure of some shale drapes (White and Willis, 2000). The estimates of the mean and the histograms (Fig. 6A) are for two-dimensional observations, and have not been corrected to three dimensions (Geehan and Underwood, 1993). The significance of the anisotropy in termination frequency can be tested statistically. The null hypothesis is that both cliff faces have the same termination frequency, based on pooling observations from both cliff faces (t = 0.37 m-1). The sampling distribution of t is computed via Monte Carlo simulation. There is greater than 10 percent

Effects of Shales in Fluvial-Deltaic Sediments

9

chance that the true t values for each of the two orientations are equal to the pooled estimate, and the observed difference is sampling fluctuation (Fig. 6B). Therefore, the null hypothesis is not rejected, and the distributions in the two orientations have not been demonstrated to be different. Horizontal anisotropy of shale dimensions is difficult to observe or quantify from two-dimensional outcrop data.

Shale Fraction and Trends In this study, shale drapes are idealized as infinitely thin but perfectly sealing barriers positioned on accretion surfaces. In this context, “shale fraction” refers to the fraction of a surface covered by shale, not the volume fraction of shale. Mudstone bodies interpreted from GPR data are modeled separately. In the two outcrop exposures, an average of 21 percent of the accretionary surfaces is covered by shale. However, the range of coverage is quite wide: about one-quarter of the surfaces have no shale, and about 20 percent of the surfaces are more than 40 percent covered by shale (Fig. 7A). Because the shale coverage is much less than 100 percent, these shales will act as baffles rather than barriers to flow. That is, they will decrease permeability and alter patterns of fluid flow, but will not completely isolate the beds bounded by the accretionary surfaces. The potential for preservation of shales depends on their position within the depositional unit. For example, tide-influenced deltas have greater preservation of shales at the base and margins of the deltas (White and Willis, 2000). In distributary channels, the preservation of shales increases upward, away from the base of the channel. This trend is observed in plots of the shale fraction versus elevation for both cliff faces (Fig. 7B). The normalized shale fraction is the shale fraction observed at a given elevation divided by the average shale fraction through the entire interval. At the top of the deposit, the shale fraction is more than 80 percent above the mean fraction, and the shale fraction at the bottom of the deposit is only about 18 percent of the mean fraction. The trend shown in Figure 7B is

z top − z f sh ( z ) = 1.82 − 1.64 f sh z top − z bot where fsh(z) is the fractional coverage of shale at the elevation of interest z,

f sh is the average shale coverage over

the entire interval, and ztop and zbot are the elevations of the top and bottom of the depositional unit. This trend is estimated by least squares subject to the constraint that



ztop

zbot

f sh ( z )dz = f sh . The trend in shale coverage may

10

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

affect flow behavior, and should be considered in constructing flow models. Methods for including this trend and its effects are discussed below. The statistics on shale lengths and fractional coverage of surfaces could be used to construct object models. However, the data are not sufficient to establish a relationship between length and elevation nor are the histograms especially robust with only 50 observations of shales. Furthermore, the histograms provide information only about the lengths of individual shales. They provide no information concerning the spatial correlation of shales, e.g., the tendency of shales to be clustered or dispersed along surfaces. Indicator geostatistics provide an alternative approach.

Semivariogram Models Indicator semivariograms are computed along each stratigraphic surface in the outcrop faces. These variograms are for shale drapes only, and are not applicable to the thicker mudstone bodies interpreted from GPR data. These variograms are used for geostatistical simulation. Because there is a trend in shale fraction (Fig. 7B), a vertically varying truncation probability is used. The indicator semivariogram is transformed to an “equivalent” Gaussian variogram to simplify the modeling process (Matheron and others, 1987). The Gaussian semivariogram is found point-by-point by finding the covariance function that satisfies

1 γ I (h; p) = p(1 − p) − 2π



arcsin CY ( h )

0

 y 2p   dθ  exp −  1 + sin θ   

where p is the cutoff proportion (p = 0.21 for this data set), yp is the corresponding inverse of the standard Gaussian distribution for a cumulative probability of p, CY is the covariance function for the standard Gaussian variable yp, and γI is the indicator semivariogram at cutoff p (Deutsch and Journel, 1998). The Gaussian correlogram CY(h) is estimated by constrained nonlinear optimization at each h. The computed Gaussian variograms and their models are shown in Fig. 8; the parameters are given in Table 2. The variograms for all surfaces are pooled to obtain these experimental variograms. Similar to the conclusions for the histograms and termination frequencies (Fig. 6), there are indications of anisotropy. There is a suggestion of a hole effect in the N-S variogram (Fig. 8). The N-S variogram indicates slightly less correlation in this direction (parallel to paleoflow), which is the same trend suggested by the histograms (Fig. 6A) and the termination frequency (Fig. 6B). The integral range for the N-S model is slightly less than for the E-W model (5.75 versus 7.60 m),

Effects of Shales in Fluvial-Deltaic Sediments

11

conforming to the other assessments of anisotropy (although these differences have not been demonstrated to be statistically significant). The possible hole effect (which is not modeled) would not affect the integral range.

Truncated Gaussian Simulation The Gaussian semivariograms (Fig. 8) are used with a sequential Gaussian simulation algorithm (Deutsch and Journel, 1998). The trend model (Fig. 7B) is used to truncate the Gaussian simulations and obtain shale indicator images:

1, y p ( x, z ) ≤ f sh ( z ) i ( x, z ) =  0, y p ( x, z ) > f sh ( z ) where yp is the cumulative probability of the standard Gaussian variable y which is simulated. The elevation for a given surface is obtained by referring to the stratigraphic framework at the location (x, y). Different trend models are considered, as discussed in the next section. Fluctuations of the mean shale fraction on the surfaces (Fig. 7A) are reproduced via fluctuations in the stochastic simulations. The indicator simulations are then used to create permeability barriers at the appropriate locations in the gridded models. Six surfaces are created for grid A, and ten for grid B (for each realization of the complete models A and B). Examples of simulated surfaces for grid A with various geostatistical parameters are shown in Figure 9.

EXPERIMENTAL DESIGN FOR STUDY OF FLOW RESPONSES Experimental design is a method to select combinations of conditions to assess the effects of process factors (Box et al., 1978). These methods have been used in reservoir engineering for sensitivity analysis (Kjønsvik and others, 1992), uncertainty assessment (Bu and Damsleth, 1996), upscaling (Narayanan, 1999), and regression (Eide and others, 1994).

Two-level Four-factor Experimental Design Four factors are selected for study (Table 3). The selected design is a full two-level factorial. Thus, 16 sets of factors (24 = 16) are considered. Two-level full factorials are simple to design and to analyze (Box et al., 1978). The factorial design allows estimation of interactions of effects. For example, the effect of shale fraction may increase as the variogram range increases. The four factors considered are (Table 3):

12

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

1.

D, the variogram range in the dip direction (aligned with paleoflow; Fig. 3). The inferred variogram ranges for all structures (Table 2) are used for the low case, and all ranges are doubled for the high case.

2.

S, the variogram range in the strike or paleoflow-perpendicular direction. It is varied in the same way as D.

3.

F, the mean fraction of the accretionary surfaces covered by shale drapes. The low value of F is the observed shale coverage fraction (0.21) and the high value of F is twice the observed shale fraction.

4.

T, the trend in shale fraction. The low value of T is 1.0, and its high value corresponds to the observed value of 1.8. The shale fraction trend is parameterized based on the maximum shale fraction, T:

 z top − z  f sh ( z ) = F T − 2(T − 1)  ......................................................(1) z top − z bot   Equation (1) results in higher fractions of shale coverage shale higher in the deposit for T > 1. The physically reasonable range is 0 ≤ T ≤ 2 . Ranges of factor variability are chosen to enclose the feasible range of the variables and to assess the effects of uncertainty in variogram parameters. All factors have approximately a two-fold range. The factor ranges allow a quantitative assessment of sensitivity for geostatistical parameters approximately equal to those observed at Corbula Gulch. The effects of varying geostatistical parameters for the shale images are shown in Figure 9. Compared to the base case (Fig. 9A), the effects of longer correlation range (Fig. 9B) and increased shale coverage (Fig. 9C) are as expected: longer correlation ranges produce larger shales, and higher shale fractions yield greater coverage of the surface. The influence of imposing a trend (Fig. 9D) is subtler, but on average the shale fraction will increase toward the left, where the accretionary surfaces are higher (Fig. 4). Because the average dips are relatively low (2.4 to 7.6 degrees), the effect of a trend in shale fraction is small for any given surface; the third surface from the bottom of grid A (Fig. 4A) is shown in Figure 9. Because different surfaces have different elevation trends (Fig. 4), the shale trends are different on these surfaces as well. For each realization of model A, six shale surfaces are created geostatistically, one for each stratigraphic surface interpreted from the GPR survey (Fig. 4). Similarly, each realization of model B requires ten shale surfaces.

Analysis Approach Because the experimental factors (Table 3) are parameters for a stochastic model for shale distributions, multiple realizations must be considered for each factor combination. In this study, 10 realizations are created for each combination of geostatistical parameters. Thus, 160 stochastic models are prepared for each stratigraphic

Effects of Shales in Fluvial-Deltaic Sediments

13

model. Grid A contains six shale surfaces and grid B contains 10 surfaces. For each combination of geostatistical parameters, the mean responses and variances among realizations is computed. These responses are then examined with response surface models and sensitivity analysis (Myers and Montgomery, 1995).

FLOW MODELING The effects of shales are assessed via numerical reservoir simulation. Reservoir simulation data sets are constructed with different images of shale distribution. Flow is simulated through these models in all three coordinate directions, aligned with the grid of the GPR surveys. The grids are oriented differently with respect to paleoflow (Fig. 3). For grid A, the x-direction is approximately perpendicular to paleoflow, and the y-direction is opposite to paleoflow. For grid B, the x- and y-directions are oblique with components along paleoflow. The zdirection is always vertical, not perpendicular to bedding. The predicted flow behavior is analyzed to quantify the importance of the parameters describing shale distribution.

Model Description Flow is simulated with a commercial reservoir simulator (Schlumberger Technology Co., 1997). The displacement process is ideal tracer flow; there are no buoyancy, capillary, relative permeability, or viscosity contrast effects (Calhoun, 1968). There are several reasons for choosing a tracer displacement rather than a waterflood as the model process. First, the tracer displacements are much less time-consuming to simulate. Second, tracer simulations isolate the effects of permeability heterogeneity and make the results less complex to interpret. Fewer factors must be considered – for example, sweep efficiency is not sensitive to rate for tracer displacements with no intrinsic dispersion. Third, truncation error can be reduced for these fully miscible systems (Rubin and Blunt, 1991). On the other hand, tracer simulations cannot investigate effects such as the interaction of gravity with reduced vertical permeability caused by shales, relative permeability and capillarity. Flow is simulated in the x-, y- and z- directions. The models are initially saturated with tracer-free fluid at constant potential (initially quiescent). For x-direction flow, the x = 0 and x = Lx faces are held at different and constant pressures throughout the simulation (Lx is the length of the grid in the x-direction). The other faces are considered impermeable, which corresponds to a no-flux boundary condition on the y- and z-faces. These boundary conditions are clear in the tracer concentration displays (Fig. 10). Fluid is injected and produced at a constant rate of 0.001 pore volumes per day. The injected fluid is marked with an ideal tracer at unit concentration. The pressures

14

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

and rates are managed via the reservoir simulator well model. For both wells, the boundary condition is infinite conductivity and nonuniform flux. Each “well” is a set of connections to gridblocks that span the entire face of the grid perpendicular to the flow direction. Because these stratigraphic grids have many layers that pinch out, there are many void blocks in these grids. Therefore, the well-to-grid connections are described explicitly rather than using automatic well-completion features in the flow simulator. This is especially important for vertical flow, where the tops and/or bottoms of the grids are defined by gridblocks in many different simulation layers. To ensure that the effects of nonhorizontal flow are included correctly in effective property calculations, the injection and production wells have the same pressure datum. The simulations for y- and z-directional flow are carried out similarly (Fig. 10). The intrinsic or “rock” x- and y-permeabilities are assumed equal, and kx / kz is assumed to be 10. Sandstone porosity averaged 19 percent. Mudstone permeability (in both drapes and bodies) is zero. Mudstone bodies are

(

)

modeled as void blocks k x = k y = k z = 0 md , whereas shale drapes are modeled as vertical transmissibility barriers (White and Barton, 1999). The results of each simulation are used to calculate three responses: upscaled permeability, breakthrough time, and sweep efficiency. The upscaled permeability is computed directly from the flow rate and the pressure drop between the wells. Reference simulations through grids with uniform permeability are used to normalize the upscaled permeability for the variable thickness in the grids. The breakthrough time is defined as the time at which the produced fluid contains more than 1 percent tracer. Sweep efficiency is the fraction of initial fluid in place that is recovered after one pore volume of injection. All simulations end after one pore volume of injection.

ANALYSIS OF RESPONSES Summary of Responses Six different models are considered to evaluate the effect of mudstone bodies and shale drapes on flow responses: 1.

AH is based on grid A; it is relatively homogeneous, with no mudstone bodies or stochastic shale drapes, and it is deterministic.

2.

AM is based on grid A and includes mudstone bodies but no stochastic shale drapes; this model is deterministic.

3.

AD is based on grid A and includes stochastic shale drapes but no mudstone bodies.

4.

AMD is based on grid A and includes mudstones and stochastic shale drapes.

Effects of Shales in Fluvial-Deltaic Sediments

5.

15

BH is based on grid B; it is relatively homogeneous, with no mudstone bodies or stochastic shale drapes, and it is deterministic.

6.

BD is based on grid B and includes stochastic shale drapes but no mudstone bodies.

As described above, 16 different combinations of geostatistical parameters (D,S,F,T) are considered for flow in the three models containing stochastic shale drapes (AD, AMD, and BD). Flow is also simulated through the models without drapes (AH, AM, and BH). The models without drapes are deterministic and have no variable factors. The drapes-only models are considered to isolate the effects of shales draping accretionary surfaces. The drapes-andbodies model (AMD) is the most realistic interpretation of the radar response of the Corbula Gulch deposits (Corbeanu, 2001). The mean responses (based on 10 realizations) for all models are presented in Table 4. The parameter combination corresponding to the observed factor values is

( D, S , F , T ) = (−1,−1,−1,+1) , where –1 indicates the

low value and +1 indicates the high value for each factor. All stochastic responses reported in Table 4 are for this factor combination. Examining the responses, the different base models (upper half of Table 4) clearly have different responses, and the responses are anisotropic. Further, the mean responses of all stochastic models are different from one another. The grids including mudstone bodies (AM and AMD) have lower upscaled permeability than the corresponding models that did not include the mudstone bodies (AH and AD). According to t-tests (Table 5), the deterministic models (no stochastic shale drapes) are drawn from a different population than the models with stochastic shale drapes. That is, shale drapes make the flow behavior of the models significantly different in most cases. The small differences in flow responses of the base and stochastic models are statistically significant because the variances among the stochastic replications are small; the coefficients of variation are always less than 0.05 and often much less. The small coefficient of variation indicates that ten realizations are probably more than adequate to predict mean response of the stochastic cases. Using ten realizations provided an accurate estimate of the variance and made the t-tests more discerning. The largest differences in responses are for z-direction flow, which is expected because of the large cross-sectional area of shale drapes in this direction. In all directions, permeability is more sensitive than any other response to the presence shale drapes. Breakthrough time is often very sensitive to heterogeneities. However, in these examples the breakthrough times are controlled by the contrast in permeability between the high-permeability layers toward the bottom of the models and the lower permeability in the accretionary bedding higher in the models (Figs. 5, 10). Consequently,

16

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

there is little variation in breakthrough time for the different models considered except for z-direction flow (Table 5, Fig. 10). Sweep efficiency also shows little variation, for similar reasons. Only very small differences in tracer front geometry can be discerned in the near-horizontal displacements without and with shale drapes (compare Fig. 10A to 10B and 10C to 10D) and the effects are moderate for vertical flow (compare Fig. 10E to 10F). On the other hand, the upscaled permeabilities are sensitive to the shale drapes in the upper part of the grids, especially in the vertical direction (Tables 4, 5). The difference in responses caused by shale drapes is never more than 2 percent for horizontal flow (Table 5). The differences in the vertical direction are more significant: from 6 to 14 percent for breakthrough time, from 1 to 4 percent for sweep efficiency, and from 13 to 23 percent for upscaled permeability. The difference in all responses is smallest for the model including massive mudstones (AM compared to AMD, Table 5) because the flow restrictions in model AM are partly redundant with the shale drapes in model AD – a flow path can be blocked only once. That is, the thin shale drapes that GPR cannot resolve are partly redundant with the resolvable mudstone bodies. However, the thin shale drapes do cause significant differences in all responses, even when mudstone bodies are present (Table 5). Thus, it is necessary to include the shale drapes using a stochastic model.

Assessing Factor Importance The contribution of each of the geostatistical factors (D, S, F, and T) are examined by analyzing the factorial design. In a two-level factorial analysis, the effects and interactions computed as simple sums and differences of the responses at the experimental points (Box et al., 1978). The main effects of the factors are the two times the average slopes of the responses with respect to each factor; interactions correspond to mixed partial derivatives. The estimates of the factors and interactions can be tested for significance with t-tests (Myers and Montgomery, 1997). If the t-value is small, the effect may be a sampling fluctuation and should not be considered significant. The factorial analysis and t-tests (at 80 percent significance level) are presented in Table 6. The normalized effect (in the rightmost column) is the effect divided by the mean response. There are several striking results in Table 6. First, the geostatistical factors have little effect on breakthrough and sweep efficiency. The only exceptions are z-direction breakthrough time of models AD and AMD. Second, the shale fraction is the only statistically significant factor. Shale fraction has a significant effect on the upscaled permeability for all models except for y-direction flow of model AD. Bearing in mind that the models without shale drapes (AH, AM, and BH) and the stochastic models including shale drapes (AD, AMD, and BD) are

Effects of Shales in Fluvial-Deltaic Sediments

17

significantly different (Table 5), these results suggest that reasonable estimates of tracer flow behavior will be obtained if the shale coverage is estimated correctly. Variogram parameters and trends in shale coverage are less important. Furthermore, there is no significant interaction between any of the factors. Finally, although the differences between the models with and without drapes are smaller when massive mudstones are included (cf. AM and AMD, Table 5), the effect of shale fraction is approximately equally significant for the mudstone and drapes-only stochastic cases (Table 6). It must be emphasized that these conclusions are valid only for the range of shale coverage, stratal geometries, flow domain dimensions, correlation ranges, and process physics that are examined in this study.

DISCUSSION Magnitude of the Effects of Shales The effects of shales are not large in this system, limited to about a 20 percent change in vertical permeability, the most sensitive response. Only the fraction of the surfaces covered by the shales is statistically significant. The lack of effect for the variogram factors is partly due to the relatively low shale coverage observed in these deposits. The interactions between fractional coverage and the other factors are not statistically significant. However, these interactions indicate that the effects of all of these factors may increase as shale fraction increases. In other words, the Corbula Gulch system is not especially sensitive to the details of the shale distribution because of the low shale coverage – with so many opportunities for “leaks,” the precise geometry of the shales is unimportant. The mean absolute value of the dip of the surfaces in the Corbula Gulch models is rather low, only about 5 degrees for both grids. This dip angle is small compared to the largest dip angles considered by Jackson and Muggeridge (2000), who predict that shales with such low dips will have relatively small effects on horizontal flow. Additionally, the fractional coverage of the shales (21 percent for the base case) is lower than the inflection point of circa 50 percent observed for geometrically simpler models (see Fig. 10 in Jackson and Muggeridge, 2000). Consequently, at Corbula Gulch the effects of the shales on horizontal properties is rather small (Table 5). If the surfaces were perfectly horizontal, there would be no effect on the horizontal properties for an ideal tracer displacement. If the dips of the surfaces were larger, the impact on horizontal permeability, breakthrough, and sweep would be increased, as would the anisotropy of the responses. The effects of the shales would increase if the accretionary surfaces were more closely spaced, as often occurs in tide-influenced settings (Willis and others, 1999).

18

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

Systems with lower intrinsic horizontal sandstone permeability, more shales, or larger shales would also be affected more strongly. The Corbula Gulch models differ from the models of Jackson and Muggeridge (2000) in several ways. In the Corbula Gulch study, three dimensional models are used throughout, sandstone permeability is nonuniform, shale placement is stochastic (using outcrop-derived geostatistics) rather than imposing regular arrays of shales, the shale dimensions and the stratigraphic frameworks are anisotropic, and the grid is stratigraphic rather than Cartesian. The Corbula models include more stratigraphic, petrophysical, and geostatistical data and therefore assess shale effects more realistically. On the other hand, the parametric models used by Jackson and Muggeridge (2000) allow analysis and generalization via dimensionless groups, which is difficult to do for the complex property distributions, stochastic features, and irregular grids used in this study. In addition, Jackson and Muggeridge (2000) considered water-displacing-oil flow, which is physically more complex than the ideal tracer flow considered in the Corbula Gulch models. Given these differences, the results of the Corbula study conform well to the predictions of Jackson and Muggeridge (2000). Shales are less significant as barriers in three dimensions than in two. Many studies of shales use twodimensional flow models (e.g., White and Barton, 1999; Jackson and Muggeridge, 2000; Willis and White, 2000); these two-dimensional models have perhaps led to overestimates of the effects of shale drapes and other heterogeneities. Further, many studies use shale “objects” rather than covariance models (e.g., Haldorsen and Lake, 1984). The objects are usually less “leaky” than covariance-based models. Desbarats (1987) uses covariance-based models in three dimensions in his study of shales and observes more significant effects than we observe in this study, but there are several differences between his models and the models considered in this study. Desbarats assigned shale properties to blocks rather than to thin interblock surfaces (as in this study), his models had more shales, and his grid is strictly Cartesian.

Upscaled Permeability and Anisotropy Although the small-scale or intrinsic permeabilities are isotropic in the horizontal plane, the upscaled x- and y-directional permeabilities are anisotropic (Table 4). The effect is most pronounced for grid A, which is aligned with paleoflow. The dip of the surfaces in grid A is anisotropic. The mean absolute dip perpendicular to paleoflow (x-direction) is 7.4 degrees, whereas the mean absolute dip parallel to paleoflow (y-direction) is 2.8 degrees. The anisotropy in dip is reflected in significant anisotropy in the flow responses (Table 4). The higher average dip occurs

Effects of Shales in Fluvial-Deltaic Sediments

19

along point-bar surfaces that accreted into the distributary channel, perpendicular to paleoflow. The paleoflow-



ky



< 1.5  : upscaled parallel (or y-direction; Fig. 3) permeability is always highest for the grid A models, 1.3 < kx   permeability is higher along accreting beds compared with across them. On the other hand, because grid B is oriented obliquely to paleoflow, the dips along the x- and y-directions are approximately equal (approximately 5 degrees), and the flow responses for models based on grid B are less anisotropic than models based on grid A

 ky   ≈ 0.9  . Recovery efficiency and sweep efficiency have anisotropy similar to that of upscaled permeability for  kx  each set of models. The upscaled horizontal anisotropy is nearly identical for models with and without shale drapes. Thus, the upscaled horizontal anisotropy is related primarily to the inclined strata, which cause variations in thickness of the high- and low-permeability regions (Figs. 4 and 5). The assumed intrinsic or rock kx / kz is 10. The upscaled horizontal-to-vertical anisotropy for the models without shale drapes is quite similar to the intrinsic anisotropy, except that the anisotropy is greater if mudstones are included. Adding shale drapes further decreases the vertical permeability (Tables 4, 5). Permeability can be estimated from analytic formulae rather than flow simulations. Cardwell and Parsons (1945) established that the harmonic and arithmetic means are the bounds on the upscaled permeability. Li and others (1999) found that tighter bounds could usually be computed. The lower bound is computed by taking the harmonic average along the flow direction for a single column of blocks, and then computing the arithmetic average over all columns. Conversely, the upper bound is computed by taking the arithmetic average for all planes perpendicular to the flow direction and then applying the harmonic average over all planes. The method can be modified for use on cornerpoint grids, subject to the same assumptions made by Li and others (1999). The arithmetic averages and estimated bounds for the models without mudstone bodies are shown in Table 4. The arithmetic average is near the upper bound except in the z-direction. In the z-direction the areal extent of permeability heterogeneities is more significant, and the estimated upper bound is less than the arithmetic average (as predicted by Li and others, 1999). The upscaled permeabilities computed by solving the steady state flow equation (in the response column of Table 4) are the most accurate. Comparing the pressure-solver results with the estimated bounds shows that the bounds do not always work well. The z-direction bounds are too low when mudstone bodies are included in the model (model AM) and too high when no mudstone bodies are included (models AH and BH).

20

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

Contrary to the results of the more accurate pressure-solver method, the bounds do not predict significant horizontal anisotropy. Fortunately, the cost of direct pressure solver calculation of upscaled permeability is not excessive for these problems.

Computational Issues Three different geologic models are considered, each with three flow directions, 16 sets of geostatistical parameters, and 10 realizations. Including additional runs for calibration, almost 1500 simulations were run to analyze the effects of the stochastic shale drapes. Such large numbers of runs present several challenges. The geostatistical models must be specified, simulated, and their parameters archived. Scripts and spreadsheets create input files for geostatistical software (Deutsch and Journel, 1997), run the software, and prepare simulator input files from the geostatistical output. All files are named based on the variable being tabulated, the grid being simulated, parameter settings, and realization number. The flow simulations are performed on approximately one-dozen desktop computers (500 MHz Pentium 3 processors). Task management software spools pending jobs to the computers when the computers are not being used. This allowed the large suite of flow simulations to be completed in only a few days. Because the runs are named systematically, they can be gathered from the networked computers and summary files are assembled on a single computer for analysis. Scripts and spreadsheets parse simulator output, compute summary statistics, and analyze the reservoir simulation results. The most time-consuming aspects of this study were the construction of the stratigraphic framework, statistical analysis of radar and fluid properties (Corbeanu, 2001), gridding, and geostatistical analysis of shale occurrence. The flow simulations were computationally intensive: the 1,464 flow simulations required 33.28 CPUdays, distributed over twelve 500 MHz Pentium 3 computers. In comparison, the efficiency of geostatistical simulation software and factorial analysis procedures is unimportant.

Future Work and Other Modeling Issues The responses might be more variable if the GPR volumes had been more areally extensive. At the scales examined, all layers except the topmost strata in grid B extend all the way across the grids (Fig. 4). If the layers did not extend all the way across the grids, then the effects of shales on the surfaces separating the strata would have impeded flow and had a greater impact on flow behavior. The intermediate-scale GPR survey acquired at Corbula Gulch (Corbeanu, 2001) could prove useful in examining these effects. Although the intermediate grid is a 2-D survey with relatively widely-spaced lines (10 m between lines, 0.25 m between traces), it covers a large enough

Effects of Shales in Fluvial-Deltaic Sediments

21

area (100 m by 100 m) to include surfaces terminating against the top and bottom bounds of the deposits. Building a flow model from this survey requires more subjective interpretation and/or stochastic modeling than is required for the smaller, denser 3-D surveys used in this study (grids A and B, 0.5 m between lines and 0.5 m between traces). Because the responses are insensitive to variogram parameters, it is reasonable to speculate that the flow responses of these deposits could be represented accurately via object models. Object models can easily reproduce the required shale fraction and that is the only significant parameter in this system. This study relies on two nonorthogonal outcrop faces to estimate the shale indicator variograms. Methods to infer or at least validate shale indicator variograms from GPR surveys would be useful. However, this task is difficult because the thickness of the shale drapes is often below the resolution of 100 MHz GPR data. Higherfrequency data, although limited in depth of penetration, could be used to attempt to obtain GPR data for shale occurrence. The chosen geostatistical methods should incorporate uncertainty in the correlation of radar properties to lithology and rock properties, perhaps using a cosimulation approach (Gómez-Hernández and Journel, 1993). The results of the factorial analysis (Table 6) indicate that shale was the only significant effect. More efficient designs can sometimes be used to reconnoiter factor effects, reducing computing requirements. A partial factorial design (24-1; Box et al., 1978) would have required only half as many simulations. However, a 24-1 partial factorial design is a Resolution III design; it confounds two-term interactions with one another. The partial factorial could not be used to determine whether any two-term interactions were significant, and determination of two-term interactions was one of the goals of this study. If more factors were being considered, a Resolution IV partial factorial design would have been an accurate, efficient choice (Box et al., 1978).

Implications for Reservoir Models The flow models for the point bar deposits at Corbula Gulch demonstrate that the effects of thin draping shales on horizontal flow are relatively minor, at least at the scale of the 3-D GPR surveys (up to 51.5 m). The shale drapes and mudstone bodies have significant effects on upscaled vertical properties, especially permeability. However, even in the vertical direction the effects of the shale drapes are moderate compared to effects reported in other investigations (e.g., Desbarats, 1987; Begg and King, 1985). The decrease in vertical permeability caused by mudstone bodies and drapes may affect gravity-influenced flow processes (Jackson and Muggeridge, 2000). Variogram and trend parameters have no significant effects on flow responses in these models. The fraction of accretionary bedding surfaces covered by shale is the only significant factor. These results can guide modelers in

22

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

constructing models for shale distribution in similar depositional settings. In analogous settings, data acquisition should focus on estimating the shale fraction with lesser emphasis on trends and ranges of correlation. Because the effects were modest in the flow models, simple models for shale distribution are likely to be adequate in analogous point bar deposits. These simplifications would not necessarily be appropriate if the dip of the accretionary surface were higher, the shale coverage were greater, shale drapes were more closely spaced, or if more of the dipping surfaces terminated against the base or top of the reservoir interval.

RESULTS AND CONCLUSIONS Shale distributions within the point-bar deposits at Corbula Gulch are described via histograms and variograms. The analysis uses data from outcrop exposures adjacent to ground penetrating radar surveys and core holes. Although a slight anisotropy is observed in all shale statistics, the anisotropy is not statistically significant. The average shale length based on termination frequency is approximately 5.4 m. Shales cover 21 percent of the accretionary surfaces at this locality. Truncated Gaussian simulation reproduces the spatial correlation of shale on the accretionary surfaces and the trend of upward-increasing shale preservation. Shales are placed on accretionary surfaces that are interpreted from the ground penetrating radar surveys. This approach is very different from most previous shale models, in which shales are placed either horizontally or with constant dip. Stochastic flow models based on ground-penetrating radar surveys and outcrop data demonstrate that for the sandstone-rich rocks at Corbula Gulch the effects of thin draping shales on flow behavior are statistically significant. Vertical flow is most affected, with shales decreasing breakthrough time (approximately 10 percent), sweep efficiency (less than 5 percent), and vertical permeability (approximately 20 percent). However, the details of the shale distribution are relatively unimportant for the Corbula Gulch deposits, at least at the scale of this study. Only the shale fraction needs to be estimated accurately.

ACKNOWLEDGMENTS This work was funded by the U.S. Department of Energy under contract DE-FG03-96ER14596 with supplementary support from BP and Chevron, corporate sponsors of the University of Texas at Dallas Sedimentology Consortium, and from The Craft and Hawkins Department of Petroleum Engineering at Louisiana State University. Adam Koesoemadinata assisted with measurements of shale dimensions from photographs of the

Effects of Shales in Fluvial-Deltaic Sediments

23

outcrop. Geologic modeling software was provided by Landmark Graphics. Flow modeling software was provided by Schlumberger Technology Co. This paper is Contribution No. __ from the Geosciences Department at the University of Texas at Dallas.

REFERENCES Armstrong, R.L., 1973, Sevier Orogenic Belt in Nevada and Utah: Geological Society of America Bulletin , v.79, 429-458. Barton, M. D., 1994, Outcrop characterization of architecture and permeability structure in fluvial-deltaic sandstones within a sequence-stratigraphic framework, Cretaceous Ferron Sandstone, Utah [unpublished PhD dissertation]: The University of Texas at Austin, 260 p. Begg, S. H., and P. R. King, 1985, Modelling the effects of shales on reservoir performance: calculation of effective permeability: 1985 Society of Petroleum Engineers Reservoir Simulation Symposium, Dallas, TX, Feb. 1013, paper SPE no. 13529. Begg, S. H., and P. R. King, 1985, Modelling the effects of shales on reservoir performance: calculation of effective vertical permeability: Society of Petroleum Engineers paper 13529, p. 331–344. Box, G. E. P., W. G. Hunter, and J. S. Hunter, 1978, Statistics for experimenters: John Wiley and Sons, New York, 653 p. Bu, T., and E. Damsleth, 1996, Error and uncertainties in reservoir performance predictions: SPE Formation Evaluation, v. 11, no. 3, p. 194-206. Cardwell, W. T., and R. L. Parsons, 1945, Average permeability of heterogeneous oil sands: Transactions of AIME, no. 160, p. 34-42. Calhoun, T. G., and R. M. Tittle, 1968, Use of radioactive isotopes in gas injection: 43rd Annual Fall Meeting of the Society of Petroleum Engineers, paper SPE 2277, Houston, TX, 29 September – 2 October. Corbeanu, R. M., 2001, Detailed internal architecture of ancient distributary channel reservoirs using groundpenetrating radar, outcrop and borehole data - case studies: Cretaceous Ferron Sandstone, Utah. [unpublished Ph.D. dissertation]: The University of Texas at Dallas, 122p. Desbarats, A. J., 1987, Numerical estimation of effective permeability in sand-shale formations: Water Resources Research, v. 23 no. 2, p. 273-286. Deutsch, C. V., 1997, Direct assessment of local accuracy and precision, in E. Y. Baafi and N. A. Schofield ,eds., Geostatistics Woolongong ’96: Kluwer Academic Publishing, Dordrecht, Netherlands, p. 115-125. Deutsch, C. V., and A. G. Journel, 1998, GSLIB Geostatistical library and User’s Guide: Oxford University Press, Oxford, United Kingdom, 369 p. Eide, A. L., L. Holden, E. Reiso, and S. I. Aanonsen, Automatic history matching by use of response surfaces and experimental design: 4th European Conference on the Mathematics of Oil Recovery, Roros, Norway, 7-10 June, 1994. Gardner, M. H., 1995, Tectonic and eustatic controls on the stratal architecture of mid-Cretaceous stratigraphic sequences, central western interior foreland basin of North America: in S. L. Dorobeck and G. M. Ross (eds.), Stratigraphic evolution of foreland basins: SEPM Special Publication no. 52, p. 243-281. Geehan, G., and Underwood, J., 1993, The use of length distributions in geological modeling, in S. S. Flint I. D. and Bryant, eds., The geological modelling of hydrocarbon reservoirs and outcrop analogues, Special Publication 15 of the International Association of Sedimentologists: Blackwell Scientific Publications, Oxford, p. 205-212.

24

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

Gómez-Hernández, J., and A. G. Journel, 1993, Joint sequential simulation of multiGaussian fields, in A. Soares, ed., Geostatistics-Troia: Kluwer, Dordrecht, Netherlands, p. 85-94. Haldorsen, Helge H; and L. W. Lake, 1984, A new approach to shale management in field scale simulation models: Society of Petroleum Engineers Journal, August, p. 447-457. Hammon, W. S., III, Zeng, X., Corbeanu, R. M., and McMechan, G. A., 2001, Estimation of the spatial variability of fluid permeability from surface and borehole GPR data and cores, with a 2-D example from the Ferron sandstone, Utah: Geophysics, submitted. Jackson, M. D., and A. H. Muggeridge, 2000, Effects of discontinuous shales on reservoir performance during horizontal waterflooding: Society of Petroleum Engineers Journal, v. 5, no. 4, p. 446-455. King, M. J., and M. Mansfield, 1999, Flow simulations of geologic models: Society of Petroleum Engineers Reservoir Engineering and Evaluation, August, v. 2, no. 4, p. 351-367. Kjønsvik, D., T. Jacobsen, and A. Jones, 1992, The effects of sedimentary heterogeneities on production from shallow marine reservoirs – what really matters?: 1992 European Petroleum Conference, London, paper SPE no. 28445. Li, D., B. Beckner and A. Kumar, 1999, A new efficient averaging technique for scaleup of multimillion-cell geologic models: 1999 SPE Annual Technical Conference and Exhibition held in Houston, Texas, 3-6 October, paper SPE no. 565554. Matheron, G., H. Beucher, C. de Fouquet, and A. Galli, 1987, Conditional simulation of the geometry of fluviodeltaic reservoirs: 62nd Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, Dallas, TX, September 27-30, paper SPE no. 16753. Myers, R. H. and D. C. Montgomery, 1995, Response surface methodology: process and product optimization using designed experiments: Wiley, New York, 700 p. Narayanan, K., 1999, Applications for response surfaces in reservoir engineering, MS Thesis, The U. of Texas at Austin 90 p. Ponting, D. K., 1992, Corner point geometry in reservoir simulation, in P. R. King, ed., Proceedings of First European Conference of the Mathematics of Oil Recovery, Cambridge, United Kingdom, p. 45-65. Rubin, B. and M. J. Blunt, 1991, Higher-order implicit flux limiting schemes for black oil simulation: 1991 Society of Petroleum Engineers Symposium on Reservoir Simulation, Anaheim, California, February 17-20, paper SPE no. 21222. Ryer, T. A., 1983, Deltaic coals of the Ferron Sandstone Member of the Mancos Shale: Predictive model for Cretaceous coal-bearing strata of the Western Interior: American Association of Petroleum Geologists Bulletin, v. 65, no. 11, p. 2323-2340. Schlumberger Technology Co., 1997, Eclipse 100 97A Reference Manual, Oxfordshire, UK, 766 p. Visser, C. A., and A. G. Chessa, 2000, A new method for estimating lengths for partially exposed features: Mathematical Geology, v. 32, no. 1, p. 109-126. White, C.D., and Barton, M. D., 1999, Translating outcrop data to flow models, with applications to the Ferron sandstone: Society of Petroleum Engineers Reservoir Evaluation and Engineering, v. 2, no. 4, p. 341-350. White, C. D., and B. J. Willis, 2000, A method to estimate length distributions from outcrop data: Mathematical Geology, v. 32, no. 4, pp. 389-419. White, C. D., B. J. Willis, K. Narayanan, and S. P. Dutton, 2000, Identifying controls on reservoir behavior using designed simulations: 2000 SPE Annual Technical Conference and Exhibition, Dallas, Texas, October 1-4, paper SPE no. 62971,. Willis B.J., J. P. Bhattacharya, S. L. Gabel, and C. D. White, 1999, Architecture of a tide-influenced delta in the Frontier Formation of central Wyoming: Sedimentology, v. 46, p. 667–688. Willis, B. J., and C. D. White, 2000, Quantitative outcrop data for flow simulation: Journal of Sedimentary Research, v. 70, p. 788–802.

Effects of Shales in Fluvial-Deltaic Sediments

25

Xu, X., 2000, Three-dimensional virtual geology: photorealistic outcrops, and their acquisition, visualization and analysis [unpublished PhD Dissertation]: The University of Texas at Dallas.

26

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

FIGURE CAPTIONS Figure 1. Section illustrating Ferron sandstone stratigraphy (modified from Gardner, 1995). The section location is shown in Figure 2. Figure 2. Location of the Corbula Gulch study site. The section X-Y is shown in Figure 1. Figure 3. Detail of Corbula Gulch study site. The location within the San Rafael Swell is given in Figure 2. The cliff faces were used to obtain data on shale distribution. The north-south (N-S) face is approximately parallel to paleoflow and the east-west (E-W) face is approximately paleoflow-perpendicular. The two three-dimensional ground penetrating radar surveys are shown as shaded areas. Additional GPR surveys and boreholes were obtained in this area but are not shown here. Points labeled CGx are locations of measured sections. The z-direction is vertically downward for both grids. Figure 4. Stratigraphic frameworks for grids A (A) and B (B). Colors indicate different stratigraphic units interpreted from GPR data (Corbeanu, 2001). Stratigraphic surfaces separate the regions with different colors. Six surfaces were interpreted for grid A, and ten for grid B. The vertical dimension shown is an average over the entire grid. Shale drapes are placed only on these surfaces. Grid location and orientation are shown in Figure 3, and related data are given in Table 1. Figure 5. Permeability models without mudstone bodies for grids AH (A) and BH (B), and with mudstone bodies for grid AM (C). The horizontal permeability ranges from 30- to 50-md, with the highest values being concentrated toward the bottom of both models. The simulation grid is constructed along the inclined strata interpreted from the ground penetrating radar survey (Fig. 4). Permeabilities were estimated using cluster analysis of GPR attributes (Corbeanu, 2001). The indicated vertical dimension is the average thickness over the entire area of the grid. The mudstone bodies have zero permeability and are considered to be inactive or “void” in the flow simulation. Mudstone bodies appear as blanked areas within grid AM (C). Figure 6. Shale length observations from outcrop faces (Fig. 3). (A) Cumulative frequency plots and Erlangian models indicate that the mean length is approximately 5 m. The observed lengths in the two cliff faces are only slightly different. (B) A Monte Carlo simulation showed that the observed difference in termination frequency, t, is not significant at the 10% level. Figure 7. Fraction of accretion surfaces covered by shale. (A) More than 25 percent of the surfaces have no shale, 20 percent of surfaces have fractional coverage greater than 0.40; 0.72 is the highest coverage fraction observed. The mean coverage is 21%. (B) The amount of coverage increases upward because of the greater potential for preservation away from the base of the distributary channels. Figure 8. Semivariogram for shale occurrence. This semivariogram is for the Gaussian variable with a cutoff corresponding to the shale fraction; it is closely related to an indicator variogram. The data are slightly anisotropic and there is a possible hole effect in the north-south (N-S) direction. The variogram pools data from all surfaces on both cliffs. Variogram parameters are given in Table 2. Figure 9. Images of shale drape occurrence for the third surface from the bottom of grid A (Fig. 4A) with varying geostatistical parameters: (A) with observed ranges, mean coverage fraction, and no trend; (B) with ranges doubled; (C) with mean coverage fraction doubled; and (D) with shale fraction increasing with elevation. Figure 10. Tracer concentration distributions at 0.5 pore volumes injected for grid A. The color scale represents tracer fraction. The dimensions of all models are those given in (A). The deterministic nodrapes model AH and a stochastic model including shale drapes AD with factors (D,S,F,T) = (-1,-1,-1,1) are shown. Flow is shown in the x- (A and B), y- (C and D) and z-directions (E and F). Except for zdirection flow, there are only minor differences in the tracer front shapes for the displacements (Tables 4, 5).

Effects of Shales in Fluvial-Deltaic Sediments

27

TABLES

Table 1. Three-dimensional flow models Model dimensions Ave. Includes X Y Grid thick. name mudstone Model Dir. of length length (m)b bodies Name +X axis (m) (m) A B

TRUE

AM

FALSE TRUE FALSE

AH BM BH

E16.0N

27.0

31.0

11.65

E23.5S

51.5

28.0

11.91

Block countsa Pore volume (m3) 2704 2924 4681 5154

Nx

Ny

Nzc

27

31

95

50

28

Active blocks 48,960

53,006 84,947 133 93,630

a

X and Y, mutually perpendicular horizontal coordinates; Z, vertical coordinate Mean thickness of model c Number of computational layers with active grid blocks b

a b

Paleoflow-parallel Paleoflow-perpendicular

Table 3. Factorial design

Factor name Dip or paleoflowparallel rangea Strike or paleoflowperpendicular rangea Shale coverage fraction Shale fraction trendb a

Symbol

Table 2. Variogram model parameters Range (m) a b N-S E-W Variance Structure None inferred Nugget Spherical 6.0 15 0.55 Exponential 30 39 0.45 Integral range (m) 5.75 7.60

Range Low

High

D

1.0

2.0

S

1.0

2.0

F

0.21

0.42

T

1.0

1.8

Multiplier for ranges in Table 2 b See Equation (1)

28

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

Direction

Table 4. Mean responses from flow models for deterministic and stochastic cases. Average Estimated bounds con Response or mean response perm. avg. perm. (md) a b Minimum Maximum k (md) τ (pv) N (pv) Grid Case BT pD (md) Model AM, mudstone bodies A AH, no bodies or drapes B

A

BH, no bodies or drapes AMD, mudstone bodies and shale drapes AD, shale drapes only

B a

BD, shale drapes only

X Y Z Base cases: X do not Y include shale Z drapes X Y Z X Stochastic Y cases: Z include X stochastic Y models for Z drapes with X (D,S,F,T) = (Y 1,-1,-1,+1) Z

0.70 0.72 0.80 0.78 0.80 0.88 0.76 0.72 0.74 0.70 0.72 0.75 0.78 0.78 0.76 0.75 0.72 0.67

0.880 0.890 0.942 0.933 0.936 0.970 0.930 0.927 0.944 0.877 0.888 0.929 0.932 0.935 0.933 0.926 0.923 0.919

23.67 35.04 2.03 35.39 45.96 5.17 41.80 38.17 5.30 23.56 34.98 1.78 35.32 45.92 4.11 41.39 37.82 4.20

40.00 3.84 43.46 3.99 37.74 3.80

39.98 39.96 1.73 43.46 43.44 40.66 37.74 37.74 37.94

Averages not computed for cases with stochastic shale drapes

Computed from steady-state pressure drop in mean uniform flow (numerical) Arithmetic average weighted by cell volume c Combination of arithmetic and harmonic averages (Li, Beckner and Kumar, 1999) b

5.53 5.69 0.10 21.89 21.94 42.46 13.74 12.05 37.71

Effects of Shales in Fluvial-Deltaic Sediments

29

Grid

A

B

Model AMD compared to AM AD compared to AH BD compared to BH

Direction

Table 5. t -tests comparing base and stochastic cases.

X Y Z X Y Z X Y Z

Risk of Type I error, α, in rejecting H0: base case drawn from same population as stochastic cases (percent) NpD τBTa k 100.00 100.00 0.26 50.44 0.05 0.14 3.06 50.07 3.06

0.00 0.16 0.85 0.05 1.49 0.08 0.00 0.03 0.00

0.00 0.12 0.02 0.04 1.07 0.02 0.00 0.01 0.00

Differences between base and stochastic responses (percent) NpD τBTa k AMD X 0.00 0.24 0.47 compared Y 0.00 0.14 0.18 to AM Z 5.92 1.36 12.91 A AD X 0.00 0.10 0.21 compared Y 2.02 0.06 0.08 to AH Z 14.18 3.94 22.91 BD X 1.68 0.42 0.97 B compared Y 0.61 0.38 0.92 to BH Z 9.45 2.66 23.11 a The resolution of 0.01 pv for τBT affects this calculation.

30

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

Model AMD

AD

BD a

b

Direction

Table 6. Responses and effects from factorial analysis of flow models.

X Y Z X Y Z X Y Z

Breakthrough time, τBT (pv) Center valuea 0.700 0.720 0.706 0.775 0.782 0.764 0.747 0.716 0.673

Shale effectb

-0.109

-0.113

Sweep Efficiency, NpD (pv) Center valuea 0.876 0.887 0.915 0.931 0.935 0.933 0.926 0.923 0.919

Shale effectb

Upscaled permeability, k (md) Center Shale a value effectb 23.49 -0.006 34.93 -0.002 1.56 -0.235 35.26 -0.003 45.88 4.11 -0.201 41.39 -0.009 37.82 -0.008 4.20 -0.205

Response estimated at the midpoint for all factors; (D,S,F,T) = (0,0,0,0)

No other effects are statistically significant. This is the normalized effect, equal to the effect divided by the response at the midpoint for all factors. Entry is omitted if not statistically significant at the 20 percent level (α = 20 percent).

Effects of Shales in Fluvial-Deltaic Sediments

31

FIGURES SW X Corbula Gulch study site

Wave- or river-dominated delta Coal Shelf margin system

Cycle 7

M I

Delta plain

Delta shoreline at time of subject point bar deposition

Blue Gate Shale

Cycle 6

J

Fe

Transgressive system

rro

nS

Cycle 5

G

an

Cycle 4 A

NE Y

Fluvial channelbelt

Sequence Boundary

ds

to

ne

0

20 km

Cycle 3 Cycle 2

C Washboard ss

Clawson ss

Tunuk Shale Figure 1. Section illustrating Ferron sandstone stratigraphy (modified from Gardner, 1995). The section location is shown in Figure 2.

Bo area of detail

ok

Y Price

Cli

ffs

6

ll we

tio ec 10

Sa n

Cr

os

ss

tch P Was a

Ferron

Emery

Ra fae lS

n(

late

F ig

au

.1

)

Utah, USA

Paleoshoreline

N

70 Corbula Gulch site

Ferron oucrop exposure 0 10 km

X

Figure 2. Location of the Corbula Gulch study site. The section X-Y is shown in Figure 1.

32

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

N

N Paleoflow rosette with 351 points

x y x

Cliff face

A CG4

Well 4

Gri

well or measured section

Grid

Well 3

y

E-W

face N-S

Origins for local coordinate systems

dB

face

CG8 CG7

100 m

Figure 3. Detail of Corbula Gulch study site. The location within the San Rafael Swell is given in Figure 2. The cliff faces were used to obtain data on shale distribution. The north-south (N-S) face is approximately parallel to paleoflow and the east-west (E-W) face is approximately paleoflow-perpendicular. The two three-dimensional ground penetrating radar surveys are shown as shaded areas. Additional GPR surveys and boreholes were obtained in this area but are not shown here. Points labeled CGx are locations of measured sections. The z-direction is vertically downward for both grids.

Effects of Shales in Fluvial-Deltaic Sediments

33

A

11.91-m

11.65-m

B

27-m

31-m

y

x

51.5-

m

28-m

z

Figure 4. Stratigraphic frameworks for grids A (A) and B (B). Colors indicate different stratigraphic units interpreted from GPR data (Corbeanu, 2001). Stratigraphic surfaces separate the regions with different colors. Six surfaces were interpreted for grid A, and ten for grid B. The vertical dimension shown is an average over the entire grid. Shale drapes are placed only on these surfaces. Grid location and orientation are shown in Figure 3, and related data are given in Table 1.

34

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

A

11.91-m

11.65-m

B

27-m

51.5-

31-m

m

28-m

C

11.65-m

y

x z

Permeability (md) 30 27-m

40

50

31-m

Figure 5. Permeability models without mudstone bodies for grids AH (A) and BH (B), and with mudstone bodies for grid AM (C). The horizontal permeability ranges from 30- to 50-md, with the highest values being concentrated toward the bottom of both models. The simulation grid is constructed along the inclined strata interpreted from the ground penetrating radar survey (Fig. 4). Permeabilities were estimated using cluster analysis of GPR attributes (Corbeanu, 2001). The indicated vertical dimension is the average thickness over the entire area of the grid. The mudstone bodies have zero permeability and are considered to be inactive or “void” in the flow simulation. Mudstone bodies appear as blanked areas within grid AM (C).

Proportion of shales shorter than

Effects of Shales in Fluvial-Deltaic Sediments

A

1

35

N-S observations

0.8 0.6 E-W observations

0.4 0.2

N-S model, t = 0.41/m E-W model, t = 0.33/m

0 0

5

10

15

20

Cumulative probability t is less than

Length (m) 1 0.8

B N-S observation

10% cutoffs

Monte Carlo simulation of distribution of t Pooled Estimate

0.6 0.4 0.2

E-W observation 0 0.2

0.3

0.4 t

0.5

0.6

Figure 6. Shale length observations from outcrop faces (Fig. 3). (A) Cumulative frequency plots and Erlangian models indicate that the mean length is approximately 5 m. The observed lengths in the two cliff faces are only slightly different. (B) A Monte Carlo simulation showed that the observed difference in termination frequency, t, is not significant at the 10% level.

Fraction with coverage less than or equal to

36

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

A

1.00 0.75 0.50 0.25

0.00 0.0

Normalized shale fraction

2

0.2 0.4 0.6 Fraction of surface covered with shale

B

0.8

Data from N-S Cliff

1.5 Model

1

0.5

0

Data from E-W Cliff

0

0.2 0.4 0.6 0.8 Fraction of thickness from top of unit

1

Figure 7. Fraction of accretion surfaces covered by shale. (A) More than 25 percent of the surfaces have no shale, 20 percent of surfaces have fractional coverage greater than 0.40; 0.72 is the highest coverage fraction observed. The mean coverage is 21%. (B) The amount of coverage increases upward because of the greater potential for preservation away from the base of the distributary channels.

Effects of Shales in Fluvial-Deltaic Sediments

37

1.6

Semivariogram

N-S model

E-W data

1.2

0.8 N-S data

0.4 E-W model 0

0

10

20

30 Lag (m)

40

50

Figure 8. Semivariogram for shale occurrence. This semivariogram is for the Gaussian variable with a cutoff corresponding to the shale fraction; it is closely related to an indicator variogram. The data are slightly anisotropic and there is a possible hole effect in the north-south (N-S) direction. The variogram pools data from all surfaces on both cliffs. Variogram parameters are given in Table 2.

38

Novakovic, White, Corbeanu, Hammon, Bhattacharya and McMechan

A

B

C

D

10

0 m

see x Fig. 3

No shale Shale

y

Figure 9. Images of shale drape occurrence for the third surface from the bottom of grid A (Fig. 4A) with varying geostatistical parameters: (A) with observed ranges, mean coverage fraction, and no trend; (B) with ranges doubled; (C) with mean coverage fraction doubled; and (D) with shale fraction increasing with elevation.

Effects of Shales in Fluvial-Deltaic Sediments

39

AH: No shale drapes (deterministic)

AD: Including shale drapes (stochastic) B Flow Direction

Flow Direction

11.65-m

X-flow

A

31-m

27-m

Y-flow

C

D

Flow Direction

Flow Direction

E

Z-flow

F

Flow Direction

Flow Direction y

x

z Tracer Concentration (fraction) 0.0

0.5

1.0

Figure 10. Tracer concentration distributions at 0.5 pore volumes injected for grid A. The color scale represents tracer fraction. The dimensions of all models are those given in (A). The deterministic nodrapes model AH and a stochastic model including shale drapes AD with factors (D,S,F,T) = (-1,-1,-1,1) are shown. Flow is shown in the x- (A and B), y- (C and D) and z-directions (E and F). Except for zdirection flow, there are only minor differences in the tracer front shapes for the displacements (Tables 4, 5).