Quantifying lateral heterogeneities in fluvio-deltaic sediments using 3-D reflection seismic data: offshore Gulf of Mexico

The quantitative characterization of depositional systems yields insight into a range of fundamental and applied geological problems. Significant contributions have been made over the last two decades in the qualitative understanding of depositional systems within a sequence stratigraphic framework [Vail et al., 1977; Galloway, 1989; Van Wagoner et al., 1990]. Seismic stratigraphy has provided a means of identifying stratigraphically significant surfaces from reflection seismic and well log data that represent the boundaries of depositional "systems tracts". There is a need to study, quantitatively, the lithologic heterogeneity within these systems tracts in order to predict their internal stratigraphic architecture.

The scaling behavior of sedimentary sequences has been a research area of considerable interest. Geologists and geophysicists have focused on the statistical description of sedimentary sequences in order to gain insight into the physical processes that guided their deposition. Observations regarding bed thickness distributions in turbidites [Rothman et al., 1994] and carbonates [Drummond and Wilkinson, 1993] suggest that their scaling characteristics reflect local sedimentation dynamics as opposed to an extrinsic forcing function. Such a forcing function might be eustatic sea-level fluctuation which Hsui et al. [1993] have suggested shows scale-invariant statistics. Modeling efforts have attempted to resolve this external versus internal forcing question [Pelletier and Turcotte , in press].

There is also a basic applied need to be able to characterize and simulate the spatial complexity of sedimentary deposits. Sedimentology based methods model geometries and spatial distributions of sand bodies within specific depositional settings [ Joseph et al., 1993; Henriquez et al., 1990; Bryant and Flint, 1992]. Geostatistical and petroleum engineering approaches that characterize the spatial distribution of rock properties such as porosity and permeability [ Haldorsen and Damsleth, 1990; Hewett, 1986; Weber, 1986; Fogg et al., 1991] are mainly concerned with the "dynamic" response of the heterogeneity field to fluid flow. Geophysical approaches attempt to relate variations in the acoustic response of sedimentary layers to geology via forward and inverse modeling techniques [Walden and Hosken, 1985; Todoeschuck et al., 1990; Bortoli et al., 1993].

Previous work on the characterization of spatial variability of rock properties has focused mainly on high resolution well log and core data [e.g. Hewett,1986; Stolum,1991]. There is significant uncertainty when well log and core data are used to infer lateral spatial variability in rock properties. Relatively few studies have focused on the statistical characterization of heterogeneity at the outcrop scale [Fisher et al., 1993; Geehan, 1993; Best et al., 1995]. One difficulty with outcrop data is that it is often difficult to assess lithologic heterogeneity in three dimensions.

In this paper we explore the possibility of utilizing 3D seismic data to statistically characterize depositional systems at the kilometer scale. The three dimensional nature of the seismic data set provides fertile ground for classical studies of depositional systems. We seek to address two fundamental questions: (i) Can seismic data be used to characterize lithologic heterogeneity at the kilometer scale? and (ii) Do different depositional systems exhibit distinct statistical behavior? We focus our statistical analyses on two depositional systems: a fluvial system and a deltaic system.

We analyze fluvial and deltaic deposits of the Pleistocene Gulf Coast. First, we establish that the seismic data is an indirect indicator of lithology in the subsurface. Then, we characterize the spatial variability of vertical and horizontal well logs, and show scale invariant behavior (2-60 meters) in these data. We analyze lateral seismic data profiles and surfaces from the two depositional systems and show that these data exhibit scale invariant behavior (100-1000 meters) and directional anisotrpy in correlation parameter, b, that is related to stratigraphic fabric. Statistical analyses of seismic data from both systems indicate distinct differences in variance, sigma², and correlation parameter, b. Finally, we compare and contrast the results obtained for the two depositional systems.

We address the problem of quantitatively characterizing lithologic heterogeneity at the kilometer scale in the Eugene Island (EI) 330 field, offshore Gulf of Mexico (Figure 1). The EI 330 field is located approximately 270 km southwest of New Orleans, Louisiana. This field is a salt withdrawal, growth faulted mini-basin [Alexander and Flemings, 1995] located in the eastern part of the Texas-Louisiana Shelf within the Gulf Coast basin. The clastic sediments of the EI 330 field are delta front facies of a Pleistocene progradational delta complex. We focus on two stratigraphic intervals: a fluvial sand (DA interval) in Block 331 and a deltaic sand (GA interval) in Block 332. In this region we have abundant wireline data and a recently shot 3D seismic survey (Figure 1). In addition, the geologic evolution of this region has been extensively characterized [ Holland, et al., 1990; Alexander and Flemings, 1995; Hart et al., in press].

Figure 1: Eugene Island 330 Field: Offshore Louisiana. The study area is located approximately 270 kilometers southwest of New Orleans, LA. Basemap shows 3-D seismic data coverage, location of wells used in this work and location of cross sections. Shaded areas marked "clinoform features" and "channel features" show areal extent of extracted seismic data used in this study.

Figure 2 is a north-south dip section that shows the location of the DA and GA sands. The DA interval is a shallow (approximately 850 m subsea depth), laterally extensive fluvial channel deposit of an incised valley fill system. Alexander and Flemings [1995] describe the DA sand as laterally continuous over tens of kilometers and interpret it to be of fluvial origin. During relative sea level fall, the region was subareally exposed (resulting in the erosional base of the sand units) and during subsequent relative sea level rise thick fluvial sands back filled these erosional unconformities. Finally, with continued sea level rise these fluvial sands were flooded by marine shales.The deeper GA interval (approximately 1500 m subsea depth) records the progradation of a shelf margin delta.

Figure 2: Generalized north-south dip section A-A' indicating location of fluvial (DA) and deltaic (GA) sands. Modified from Alexander and Flemings [1995].

Hart et al. [in press] interpret the GA interval as part of a shelf margin lowstand delta complex with vertical and lateral interfingering of erosional channel, prograding clinoform and base of slope failure deposits over the kilometer scale. Channel features in the DA and clinoform features in the GA intervals are observed on seismic sections in the study area (Figures 3(A) and 3(B)).

Figure 3: Vertical cross sections of seismic data showing channel and clinoform features. (A) East-west seismic section B-B' showing interpreted channel features within fluvial (DA) horizon. Location of line B-B' indicated on Figure 1. Figure 3: Vertical cross sections of seismic data showing channel and clinoform features.(B) North-south seismic section C-C' showing clinoform features mapped within deltaic (GA) horizon. Location of line C-C' indicated on Figure 1. Progradation direction for clinoforms is from north to south.

We seek to quantify lateral lithologic heterogeneity using stacked and processed multichannel 3D seismic data. To do so, we need to establish that the seismic response is driven by variations in lithology. We address this issue by establishing the correlation between lithology and acoustic impedance from vertical well log data, and by demonstrating a correlation between synthetic seismograms (generated from acoustic impedance profiles) and actual seismic data for specific well locations. We show that seismic data in water saturated zones record changes in lithology and as a consequence can be used to characterize lateral lithologic heterogeneities at the 100-1000 meter scale.

We use a stacked, processed multichannel 3D seismic data set . The seismic survey was stacked and deterministically deconvolved to remove the source wavelet and then migrated in time in both the inline and crossline directions. The acquisition bin spacing is 22 meters inline by 15 meters crossline. A running summation (integration) has been performed trace by trace on the 3D data cube to yield a band limited acoustic impedance function. Peterson et al. [1955] show that an integration of reflection coefficients, RFC, may be approximated by one half the natural logarithm of the acoustic impedance:

where rho is the bulk density in Kg. m^-3; V is the bulk velocity in ms^-1 , rho*V is the acoustic impedance in Kgm^-2s, and z is the depth in meters. The seismic data used in this study is thus a band limited acoustic impedance function. The term "impedance function" will be used throughout the rest of this paper as defined by equation (2). The impedance function may be evaluated by integrating the reflection coefficients, RFC (calculated from log data with equation 1) or by substituting wireline-derived sonic velocity and density data into equation (2). Inferring lithological variations from actual seismic data involves relating the seismic wavefield to the acoustic impedance structure and ultimately to the geology. We demonstrate below that there is a correlation between lithology and acoustic impedance in wireline data, and between acoustic impedance function derived from sonic and density well logs and actual seismic data.

We investigate the correlation between lithology and acoustic impedance well log data in the vertical direction for the DA and GA intervals. The gamma ray (GR) log is used as an indicator of lithology and the acoustic impedance is computed using the sonic (DT) and bulk density (RHOB) logs. A cross-plot of GR versus acoustic impedance logs from well 331_SH_1 for a sand shale interval in the vicinity of the fluvial DA sand shows that, although there is an overlap of gamma ray log ranges for the sand shale sequence in the moderate acoustic impedance range (3.6-4.3 x 10 ⁶ Kgm^-2s^-1), clean sands and shales, in general, have different impedance ranges (Figure 4(A)). A cross plot of GR versus acoustic impedance from well 331_SH_2 shows similar behavior for the GA sand (Figure 4(B)). For both the log and seismic analyses we chose regions that are water saturated and thus do not consider the effect of hydrocarbons. The crossplots indicate that clean, water saturated sands have, in general, a lower acoustic impedance than the shales in this study area as is often encountered in normally compacted, relatively shallow strata in the Gulf of Mexico. Thus, in a given sand shale sequence we expect gross sand-shale facies changes to be captured by the impedance function.

Figure 4: Crossplot of gamma ray (GR) and acoustic impedance well log data for fluvial and deltaic sands. (A) fluvial (DA) interval. The data points in the crossplot are from the log interval shown at the left. Shaded area on logs indicates location of DA interval. DT is the sonic log in m-seconds/meter and RHOB is the bulk density log in Kgm-3. Data obtained from well 331_SH_1. Figure 4: Crossplot of gamma ray (GR) and acoustic impedance well log data for fluvial and deltaic sands. (B) deltaic (GA) interval. The data points in the cross plot are from the log interval shown at the left. Shaded area on logs indicates location of GA interval. Units for DT and RHOB logs is the same as in (A). Data obtained from well 331_SH_2.

Secondly, we study the correlation between the actual seismic trace and a synthetically generated impedance function derived from well log data. The synthetic well log derived impedance function is filtered to the seismic bandwidth (10-50 Hertz) to facilitate a comparison between the traces. The synthetic impedance function for well 331_SH_1 was obtained by integrating the reflection coefficients. We compare only the relative positive and negative excursions of the synthetic impedance function and actual seismic trace. The absolute values of the actual seismic trace are different from the synthetic impedance function due to operations performed in acquisition and processing of seismic data. The results indicate a good match at the abrupt sand shale interfaces for this well (Figure 5(A)). Figure 5(B) shows the correlation between actual and synthetic seismic traces in the vicinity of the GA interval for well 331_SH_1 . Summarizing the results, we have established that clean water saturated sands have lower impedances than shales for the DA and GA sections, and synthetic impedance function traces derived from well logs correlate relatively well with the actual seismic traces extracted from the 3D data cube.

Figure 5: Correlation between synthetic well log derived impedance function with actual seismic data tracefor well 331_SH_1. (A) fluvial (DA) sand over a 170 meter vertical section (0.8 to 0.98 sec onds two-way travel time). The impedance function was calculated by integrating the reflection coefficients (RFC) obtained from Equation (1) and filtered using a band pass filter to retain frequencies comparable to that present in the actual seismic trace. The logs are sampled in 1 milli-second(ms) two-way travel time. The checkshot data from well 331_SH_1 was used to construct the time-depth tie. (B) deltaic (GA) sand over a 200 meter vertical section (1.32 to 1.5 seconds two-way travel time).

Next, we characterize the spatial variability of lithologic, impedance function and porosity data from vertical and horizontal well log data. We then characterize the one and two dimensional lateral spatial variability in our seismic dataset by analyzing profiles and surfaces from within specific stratigraphic levels. Our inverse approach which assumes that stacked multi-channel seismic data can be used to infer the acoustic impedance, and ultimately lithologic, heterogeneity compliments the forward-modeling approach of Hurich [1996] who showed a strong correlation between the statistical properties of an acoustic impedance field and the resulting synthetic reflection seismic wavefield.

We use Fourier spectral analysis and variogram methods to characterize the variability in our well log and seismic data. The use of both techniques provides a more robust means of assessing the scaling behavior in our data. We outline each of these techniques below:

Spectral Analysis:
Random functions can be characterized by variogram based methods, covariance functions or Fourier spectral analysis. Spectral analysis seeks to describe the frequency content and scaling behavior of a signal. The relationship between spectral density and frequency can be used to define the scaling characteristic of any random "noisy" signal [ Voss, 1988]. Since we are dealing with random functions (such as reflection amplitude, well log response, etc) that vary with space coordinates, frequency is replaced by wave number. The spatial variability of a random function exhibits scale invariance if the spectral density S(k) has a power law dependance on wave number (k) as:

The exponent b is a measure of the variability present in a signal and is zero for white noise. The larger this exponent, the smoother and more correlated is a given signal (Figure 6).

Figure 6: Examples of one dimensional profiles in the space and wave number domains with varying degrees of correlation. b value captures the degree of "noise" in the profiles. S(f) is the power spectra as a function of frequency or wavenumber. Figure modified after Voss (1988).

If the underlying distribution is assumed or known to be Gaussian, as in the case of fractional Gaussian noise (fGn) and fractional Brownian motion (fBm) models, then the signal can be statistically defined by three parameters: the exponent,b, the mean (m) and variance (sigma²) of the distribution. For fGn, b lies between -1 and 1 and for fBm b lies between 1 and 3. Random signals with power law spectra but non-Gaussian probability density functions have been studied by Lovejoy and Schertzer [1986], Gray et al. [1993] and Painter and Paterson [1994]. The parameter b, obtained from spectral analysis of one dimensional data, is related to the fractal dimension, D, by the following relation [Turcotte, 1992]:

for a fBm model with Gaussian increments and 1<b<3. Fourier spectral methods and scale-invariant statistics have been previously used in the characterization of earth topography [ Brown and Scholz, 1985; Huang and Turcotte, 1989; Rapp, 1989; Huang and Turcotte, 1990; England, 1992], bathymetry [ Gilbert and Malinverno, 1988; Malinverno, 1989] and reservoir heterogeneity [ Hewett, 1986; Stolum, 1991, Tubman and Crane, 1995].

Variogram Analysis:
The variogram is a measure of spatial continuity that measures the correlation between pairs of data values, z(x), separated by a specified distance or lag (u). The traditional semi-variogram g(u) is defined by [Isaaks and Srivastava, 1989; Deutsch and Journel, 1991]:

where N(u) is the number of data pairs separated by lag u; z(x+u) is the data value located at position (x+u); z(x) is the data value located at position (x) The variogram provides a variance estimate as a function of spacing between data points. Traditional variogram models assume that closely spaced data values (small lags) have low variance while values separated by large lags are uncorrelated. Note that scale invariant data exhibit long range correlations and the variogram does not reach a limiting variance at larger lags. The variogram and power spectrum are related via the covariance function [Hewett, 1986]. For a fBm process a necessary, but not sufficient, relationship between the slope, s, of the log(g(u)) versus log(u) plot and the parameter, b, is b=s+1. The appendix provides a detailed description of the statistical methodology used in this study.

Figure 7: Well log data used in statistical analysis (A) Actual log traces. Data sampled every 0.3048 meters (1 foot) for well 331_SH_1 and every 0.1524 meters (0.5 foot) for wells 330_SH_B-14 and 330_PZ_D-11. Impedance function defined by equation 2. Depths are true vertical depths except for horizontal well 330_PZ_D-11 where depth refers to measured depth. Figure 7: Well log data used in statistical analysis(B) Histograms of log data.

The lithology data from the vertical wells are bi-modally distributed reflecting gamma ray and spontaneous potential values concentrated in the sand (gamma ray ~ 30 API units; spontaneous potential ~ -26 mV) and shale (gamma ray ~ 70 API units; spontaneous potential ~ -19 mV). Both the lithology and impedance function data exhibit scale invariant behavior over the 2-60 meter scale (Figure 8(A)). Interestingly, both the lithology and impedance function data show similar correlation with an average b approximately equal to 1.8 (specific values are shown in Figure 8(A) and Table 1).

We contrast observations of vertical well log statistics with that of one horizontal well (330_PZ_D-11). This horizontal well was drilled entirely within the fluvial (DA) sand. Only a lithology (gamma ray) log was available for this well and it exhibits much lower variance in gamma ray values than that observed in the vertical well 331_SH_B-14 (Figure 7(B); Table 1). The spectra of the horizontal well gamma ray log indicates a lower b value (Figure 8(A)) than that obtained for lithology logs from the two vertical wells (b = 1.54 versus b = 1.75). Our analysis of the few available gamma ray logs from horizontal wells indicates that tool or measurement noise may contribute to the lower b values. Although this indicates a noisier profile in the horizontal direction, it is the difference in variance, sigma² (269.94 versus 22.65) between vertical and horizontal lithology data that characterizes the obvious lithology anisotropy observed in these sedimentary deposits.

Figure 8: (A) Power spectra of log data. Figure 8: (B) Variograms of log data. A Hanning windowed periodogram estimator was used to compute power spectra of well log data. The variograms were computed for lags up to one-third of the data length [Hardy and Beier, 1994] since the results may not be statistically meaningful for longer lags. Slopes of the power spectra-wavenumber (b) and variogram-lag (s) plots are indicated on the figures. Slopes were determined with a least squares fit of the data plotted on logarithmic coordinates.

Next, we analyze the spectra of sonic, as well as sonic log-derived porosity and velocity, density and acoustic impedance logs from well 331_SH_1 over the same vertical stratigraphic section as the previous data. Comparing spectra of sonic, velocity and porosity data, we observe that all three have the same b value of 1.36 (Figure 9).

Figure 9: Power spectra of sonic-derived porosity, sonic, bulk density, velocity and acoustic impedance logs from well 331_SH_1. Porosity values for well 331_SH_1 were computed from the sonic log using the empirical relationship f= 0.67*(Dtlog-Dtmatrix)/Dtlog [Schlumberger, 1987].

Spectral analysis of the bulk density log reveals a smoother profile (b = 1.66) and the b value for the acoustic impedance, the product of density and velocity, lies between 1.36 and 1.66 (b = 1.57). Also, as expected, linear transformation of the sonic log data to velocity and porosity has no effect on the slope of the power spectral density-wavenumber plot. We note that spectra of density-derived porosity, obtained from linear transformation of the bulk density log, will have the same b value as the density log spectra (i.e. b=1.66). This difference in sonic and density-derived porosity spectra may be related to the effects of adverse hole conditions on tool response. However, spectral analysis of sonic log-derived porosity data over a larger (1200 meter) vertical sand shale section from well 331_SH_1 provided a b value of 1.15 indicating that the b parameter is affected by the length of the section analyzed. This effect of data length on the b value has been noted by other researchers [Hardy and Beier, 1994; Gray et al., 1993]. We also analyzed neutron porosity logs from two wells (315_EX_1 and 337_PZ_3) over a 1200 meter interval and obtained an average b value of 1.25. Pelletier and Turcotte [in press] report similar values for neutron porosity logs analyzed over a comparable vertical stratigraphic interval in the Eugene Island 330 field area (average b value=1.41).

We now place our results in the context of previous work on scaling behavior in well log data. Walden and Hosken [1985] analyzed the power spectra of reflection coefficients (as given by equation 1) from eight well logs and reported a decay in power spectrum according to a power law f^-b,, where f is frequency and b lies between -0.5 and -1.5. They show that integrating the reflection coefficients will result in decay of the power spectrum according to f^-b where b is between 0.5 and 1.5 for the impedance function which is lower than the average value of 1.8 that we obtain. Todoeschuck et al. [1990] also report scale invariant behavior in reflection coefficient data obtained from 18 wells in the Peace River district of Alberta with an average b value of -0.68 (i.e. a b of 1.32 for the integrated reflection coefficients). Holliger [1996] reported b values between 1 and 1.2 from the spectral analysis of 10 P-wave sonic logs from upper crustal drill sites in Europe and North America. Tubman and Crane [1995] and Hardy and Beier [1994] report lower b values (0.5-0.9) for neutron porosity and density logs from vertical wells. As pointed out by Pelletier and Turcotte [in press] the results obtained by Tubman and Crane [1995] are different since they indirectly computed b from the rescaled range statistical technique instead of the spectral analysis method. Hardy and Beier's [1994] lower b values may reflect a significantly different geologic setting (shallow water carbonates).

None of our b values, obtained from spectral analysis of vertical well log data, are less than 1. Vertical well log variability has been mostly modeled as fractional Gaussian noise (fGn) with an underlying Gaussian distribution and b values less than 1 [Hewett, 1986; Todoeschuck, et al., 1990; Hardy and Beier, 1994; Molz and Boman, 1993; Tubman and Crane, 1995]. Recent work by Holliger [1996], Pelletier and Turcotte [in press] and Pelletier [pers. comm.] suggest that spectral analysis of vertical log data commonly yield b values between 1 and 2 and are not necessarily restricted to being less than 1. If we assume an underlying Gaussian distribution, then this range for b values ( i.e. b&rt1) suggests that variability in our log data could be modeled as fractional Brownian motion [Mandelbrot and Van Ness, 1968; Saupe, 1991; Hardy and Beier, 1994]. A comparison of the slopes of the power spectrum (b) and variogram (s) plots (Figure 8; Table 1) shows that most logs satisfy the relationship b=s+1. As noted previously, although this relationship does not conclusively indicate fBm behavior, it does provide additional information that is in agreement with the characteristics of fBm. Gray et al. [1993] compared the value of the Hurst parameter, H, obtained from spectral analysis of Gaussian and non-Gaussian data with similar variability and showed that the H values were not considerably affected (within approximately +/- 0.1 ) by different distributions. For fBm, the relationship between b and H is b=2H+1 and for fGn it is b=2H-1 [Hewett, 1986] indicating that the scaling exponent, b, is relatively robust to changes in the underlying distribution.

We used the Capon estimator [Capon, 1969; Kay, 1988] instead of the traditional periodogram spectral estimator to compute the power spectra for the seismic data sequences primarily because of the short length of the seismic data sequences. It has been recognized that the periodogram spectral estimator does not work well with short data records and results in a spectral estimate that has a high variance [Austin et al., 1994]. This higher variance in the spectral estimate results in a greater uncertainty regarding the least squares estimate of the slope, b, of a power spectral density-wavenumber plot. The Capon estimator, a minimum variance estimator, provides a spectral estimate with lower variance and consequently a more accurate b value for short data records. For calibration, we compared spectral estimates of well log data (of sufficient data length) obtained from the Capon and modified periodogram estimators, and the resulting scaling relations were identical.

We first map the top of the DA and GA horizons on the seismic sections and then extract horizontal, two dimensional seismic data surfaces from within each interval. First, we study the degree of correlation present in one dimensional profiles and then analyze the two dimensional statistical variation of the seismic data surfaces for both depositional systems.

Results from statistical anlaysis of 1D seismic data profiles

1. Fluvial system (DA)
We map channel features within the DA interval on the basis of seismic and wireline character. The top of the DA sand forms a prominent reflector in the 3D seismic data. We extracted and displayed in map view the amplitudes on a horizon 20 milliseconds below this mapped horizon, which is roughly in the middle of the sand (Figure 10(A)). The areal extent over which we map these features is shown in Figure 1. The amplitude map reveals linear, high amplitude events trending north-northwest. These events correspond to fluvial channels and are delineated with arrows (Figure 10(A)).

Figure 10: Map view of two dimensional seismic data surfaces. (A) fluvial (DA) interval. Arrows indicate location of channel features and dashed lines show position of one dimensional profiles analyzed in this study. (B) deltaic (GA) interval. Solid white lines indicate location of clinoform features and dashed lines show position of one dimensional profiles analyzed in this study. Figure 3(B) shows the same clinoform features in vertical cross section. Dark colors are associated with negative seismic amplitudes. Note that the inverted seismic dataset used in this study provides a first order approximation of acoustic impedance variations.

The seismic data profiles are smoother than observed in the well log data for the obvious reason that these data do not have a high frequency component (Figure 11(A)). Even by visual examination, the east-west profile appears to vary more than the north-south profiles (Figure 11(A)). Histograms of these data (Figure 11(B)) show that the amplitudes are distributed around a zero mean.

Figure 11: One dimensional seismic data profiles from fluvial (DA) interval. (A) Actual seismic traces. Two traces oriented north-south (N-S) and one trace oriented east-west (E-W)as shown in Figure 10(A). North-south traces aligned along dominant channeling direction. (B) Histograms of seismic data profiles. Variance of data profiles listed in Table 2.

Spectral analysis of these one dimensional profiles (Figure 12(A)) reveals a greater correlation in the north-south profiles ( ) than there is in the east-west profile ().

Figure 12: Spectral and variogram statistics for one dimensional seismic profiles from fluvial (DA) system (A) Power spectra (B) Variograms. Capon estimator was used to compute power spectra of profiles. The variograms were computed for lags up to one-third of the data length and the slope determined from the linear portion of the plot to a lag of 400 meters. Slopes of power spectra (b) and variogram (s) plots were calculated from a least squares linear fit to the log-log plot of power spectra (variogram) versus wavenumber (lag) and are listed in Table 2.

The variance in the north-south profiles is slightly less than that in the east-west profiles (Table 2). The directional anisotropy present appears to be associated with the stratigraphic fabric; channel parallel profiles being more correlated than those that are perpendicular to the dominant channel direction.

Spectral analysis of the profiles reveals b values greater than 1(Figure 12(A)). If the underlying distribution were Gaussian, then this relationship suggests fBm behavior. A comparison of the slopes of the variogram and power spectra (Table 2) show that two of the three profiles approximately satisfy the relationship b=s+1 in agreement with the characteristics of fBm.

2. Deltaic system (GA)
We next analyze one dimensional profiles from a horizon parallel to the GA sand. Again, we mapped on the prominent seismic reflector at the top of the GA sand and extracted a seismic horizon from 30 milliseconds beneath the mapped horizon (Figure 10(B)).

Figure 10: Map view of two dimensional seismic data surfaces. (B) deltaic (GA) interval. Solid white lines indicate location of clinoform features and dashed lines show position of one dimensional profiles analyzed in this study. Figure 3(B) shows the same clinoform features in vertical cross section. Dark colors are associated with negative seismic amplitudes. Note that the inverted seismic dataset used in this study provides a first order approximation of acoustic impedance variations.

A map view of the areal extent of this extracted seismic horizon is shown in Figure 1. The amplitude map reveals three prominent, curved lineations oriented approximately east-north-east. These lineations are strong negative seismic amplitudes that correspond to prograding, sand rich clinoforms when viewed in cross section (Figure 3(B)). The extracted one dimensional profile data and associated histograms are shown in Figure 13. Both the profiles approximately parallel to strike (E-W) and those perpendicular to strike (N-S) have variances that are an order of magnitude greater than that in the fluvial (DA) profiles (Table 2 and Table 3). This reflects the fact that these delta front deposits contain interbedded sands and shales while those in the fluvial interval were sand dominated.

Figure 13: One dimensional seismic profiles from deltaic (GA) interval. (A) Actual seismic traces. Traces oriented perpendicular and approximately parallel to clinoform features shown in Figure 10(B). Figure 13: One dimensional seismic profiles from deltaic (GA) interval. (B) Histograms of seismic amplitude profiles. Note change in scale of seismic data values on histogram as compared to Figure 11(B) indicating higher variance in deltaic interval. Variance of profiles listed in Table 3.

Spectral analysis of the four one dimensional profiles reveals mild directional anisotropy (Figure 14(A), Table 3). The two east-west profiles have lower b values () than the north-south profiles (). It is somewhat surprising that there is less correlation parallel to the major clinoforms. This may be because the extracted profiles are not taken exactly parallel to the clinoforms or because the clinoforms themselves are a small proportion of the whole data set. Regardless, spectral and variogram analyses (Figure 14) reveal scale invariant behavior over the 100-1000 meter scale.

Figure 14: Spectral and variogram statistics for one dimensional seismic profiles from deltaic (GA) system (A) Power spectra Figure 14: Spectral and variogram statistics for one dimensional seismic profiles from deltaic (GA) system (B) Variograms. Capon estimator was used to compute power spectra of profiles. The variograms were computed for lags up to one-third of the data length. Slopes of power spectra (b) and variogram (s) plots were calculated from a least squares linear fit to the log-log plot of power spectra (variogram) versus wavenumber (lag) and are listed in Table 3.

As in the fluvial system, profiles from the deltaic system have b values greater than 1 (Figure 14(A)) and a comparison of the variogram and power spectrum slopes shows that most profiles approximately satisfy b=s+1 in agreement with the characteristics of fBm.

Results from the statistical analysis of 2D seismic data surfaces Spectral analysis can be extended to statistically characterize two dimensional surfaces. A first order estimate of the roughness or variation of a surface may be obtained by computing an averaged spectral estimate over a number of one dimensional profiles in the direction of interest. Figure 15(A) shows the averaged spectral parameters for seismic data from both the fluvial (DA) and deltaic (GA) depositional systems. Also, the b value in the north-south direction for the fluvial system is lower than that for the two one-dimensional profiles analyzed earlier since it represents an average value over a number of one dimensional profiles both within and outside the dominant channeling direction. The averaged one dimensional b values for the deltaic system do not indicate significant directional anisotropy as shown in Figure 15(A). A more rigorous analysis may be performed to determine the 2D spectral parameter [ Turcotte, 1992] for the seismic data surfaces from the fluvial and deltaic systems. This analysis characterizes the variation or roughness across the entire two dimensional surface and is not particularly helpful in focusing on specific features that contribute to directional anisotropy in the b value. Nonetheless, we note that the fluvial and deltaic systems exhibit distinct two dimensional b values as shown in Figure 15(B).

Figure 15: Two dimensional power spectra of seismic data from fluvial (DA) and deltaic (GA) intervals. (A) Average of one dimensional power spectral estimates. Average taken over every tenth trace extracted from the map view shown in Figure 10(A) for fluvial system and from Figure 10(B) for deltaic system. The Capon estimator was used to determine the power spectra.
Figure 15: Two dimensional power spectra of seismic data from fluvial (DA) and deltaic (GA) intervals. (B) Two dimensional spectra of seismic data surfaces shown in Figure 10(A) and 10(B) for fluvial (DA) and deltaic (GA) intervals respectively. A 2D spectral estimator [Turcotte, 1992] was used to determine the power spectra of the two dimensional surface (see Appendix).

We have shown that seismic data at the kilometer scale can be related to lateral lithologic heterogeneity within specific depositional systems. Laterally extensive channel deposits of the incised valley-fill depositional system have a higher correlation structure along the depositional dip than orthogonal to dip, revealing a lithologic anisotropy related to the geometry of the fluvial channels. In contrast, progradational deltaic deposits show approximately similar spatial correlation both along strike and dip orientations. In contrast to the fluvial deposits, these deltaics have much larger variances reflecting the interbedding of mud and sand in their deposits.

Our results have implications in a range of spheres. Hewett [1986] proposed that vertical porosity variations were characteristic of fractional Gaussian noise (fGn) and that horizontal heterogeneity was characteristic of fractional Brownian motion (fBm). Since that time there has been a range of work that has at least tacitly confirmed this proposal as far as vertical porosity data is concerned [Molz and Boman, 1993; Hardy and Beier, 1994; Tubman and Crane, 1995] yet there has been remarkably little detailed exploration of this behavior for other well log data. In contrast, spectral analysis of our well log data (gamma ray, density, sonic and porosity) indicates scaling behavior characteristic of fractional Brownian motion (fBm). Pelletier and Turcotte [in press] and Holliger [1996] also report scaling behavior characteristic of fBm (i.e. b greater than 1) from spectral analysis of porosity and sonic log data respectively. While we observe directional anisotropy in b value, both vertically and horizontally, it is mainly the variance of the underlying distribution (eg. greater in the vertical than horizontal) that contributes to the characteristic layering in Eugene Island strata.

Conventional modeling efforts which asssume an uncorrelated structure at the larger scale (such as the variogram technique) will not reproduce the long range correlations (100-1000 meter scale) we observe. Scaling models which do consider long range correlation have been extensively used to simulate subsurface variations in rock properties [Hewett, 1986; Molz and Boman, 1993; Painter, 1996]. The scaling parameter, b, is an important variable in the modeling procedure and our observations regarding the b values in fluvial and deltaic depositional settings could provide a constraint on modeling lateral heterogeneities of rock properties.

This work may have application in geologic modeling of lithologic heterogeneity within a sequence stratigraphic framework. Sequence stratigraphic analysis has evolved from defining major bounding surfaces or sequence boundaries, to describing the stratigraphic architecture (depositional systems tracts) within sequences [ Vail et al., 1977; Jervey, 1988; Posamentier et al., 1988; Van Wagoner et al., 1990]. Our work has the potential to quantitatively characterize lithologic heterogeneity at the systems-tract scale. A knowledge of this heterogeneity may assist researchers studying the connectivity of individual sand bodies within these systems and their effect on fluid transport [ Weber et al., 1995].

While we have statistically characterized the correlation structure in seismic and well log data within specific depositonal systems, we have not addressed what physical mechanism is responsible for the observed scale invariant behavior of sedimentary rock properties. A number of researchers have used probabilistic sedimentation models to predict bed thickness distributions [e.g. Kolmogorov, 1951; Dacey, 1979; Strauss and Sadler, 1989]. Pelletier and Turcotte [in review] used a stochastic diffusion model of sediment transport and concluded that successive synthetic bed thicknesses generated by their model were uncorrelated. Plotnick and Prestegaard [1995] report long range correlations in bed thicknesses from analysis of actual data of braided and meandering river deposits. If the spectral statistics of lithology indicators (such as the gamma ray log) and bed thickness are similar, then our results support the existence of these long range correlations. Continued characterization of lithologic heterogeneity in a variety of depositional settings will provide an important constraint on these modeling efforts.

Finally, as Austin et al. [1994] pointed out, power law spectra introduce certain difficulties in the spectral estimation process especially when data records are short. Our results regarding the b parameter and spectral estimates of short seismic data sequences using the Capon estimator [Capon, 1969; Kay, 1988] should be more robust than estimates obtained using conventional Fourier based methods. We have found that it is particularly important to develop a rigorous and repeatable procedure to the spectral estimation process which includes computation of the slope of the log(power spectra)-log(wavenumber) plot within a well established range of wavenumbers (see Appendix). We believe that the results of previous work related to computation of spectral slopes would be improved upon if the assumptions behind determining the slopes were clearly outlined and that this would also facilitate comparison between scaling exponents, such as b, obtained from a variety of datasets.

Well log and 3D reflection seismic data in our study area exhibit scale-invariant statistics over the meter to kilometer scale. Spectral and variogram analysis of our well log and seismic data suggests scaling behavior (b&rt1) characteristic of fractional Brownian motions (fBm).

We document directional anisotropy in the correlation parameter, b, from spectral analysis of seismic data profiles from within fluvial and deltaic depositional environments. Statistical analyses of seismic data reveals that fluvial and deltaic systems exhibit distinct statistical behavior. Our use of Capon's spectral estimator should provide more robust results regarding the b parameter than will conventional Fourier based techniques, particularly in the spectral analysis of short data sequences. A consistent approach to determining the slope of the power spectra-wavenumber plot that includes defining the lower and upper wavenumber bounds prior to line fitting will ensure repeatability of the results.

Modified Periodogram Spectral Estimator:
One straight forward technique in estimating the power spectrum estimate of a finite data series is to compute the DFT of the series and take the magnitude squared of the result. This simple estimate is called the periodogram. The problem with the periodogram is that the variance in the estimated spectrum is large and consequently this method does not provide a stable estimate of the power spectrum. Welch [1967] suggested the use of a windowing function prior to computing the periodogram in order to reduce spectral leakage effects and diminish the variance in the spectral estimate. We employ Welch's formulation in computing the power spectra of a finite data series.

Consider a one dimensional, single valued function, Z(x_n), of depth or location, x, such as a gamma ray or porosity log. Suppose we have data values measured at N locations xi (i=0,1,...N-1) spaced at intervals delta. The modified periodogram spectral estimate, , based on Welch's formulation is given by [ Kay, 1988]:

where the Hanning window is

and the normalization factor, U, is

Prior to taking the Fourier transform, the data series was detrended to remove any linear trend present in the data. The modified periodogram estimator works well for data record lengths greater than a few hundred values. The well log data records analyzed in this study had a few thousand data values and could thus be analyzed using this method.

Capon Spectral Estimator:
As demonstrated by Austin et al. [1994] traditional Fourier based methods have higher variance in the spectral estimate and this is particularly problematic for short data records. The one dimensional seismic data profiles that we analyze contain between fifty and two hundred data values. Spectral analysis of these data requires a more accurate spectral estimator that works well with short data records. We employ Capon's minimum variance spectral estimator [ Capon, 1969; Kay, 1988; Austin et al., 1994] to estimate the power spectra of the seismic data profiles.

Consider the data series Z(x_n) measured at N locations. Following the method outlined by Kay [1988] we first compute an estimate of the elements of the autocorrelation matrix, R_zz:

where p (~0.3N) is the dimension of the autocorrelation matrix. The Capon spectral estimate, P_CAP^(k), is given by

and the vector, e, is

The operators H and T are the Hermitian (or conjugate) and normal transpose operators respectively.

Estimation of Two Dimensional Spectra:
The two spectral methods discussed previously may be used to extend the analysis to two dimensional data surfaces. A simple way of estimating an "average" power spectra for a specific direction on a two dimensional surface is to average a number of one dimensional spectral estimates from that particular direction. We used this technique to estimate an "average" power spectra and b parameter in the north-south and east-west directions for our seismic data surfaces. We computed the spectral estimate of a number of one dimensional profiles using the Capon method and averaged these estimates over all the profiles. Although this method does not give the actual two dimensional spectra of the data surface, it does provide some insight into the "average" statistical variation of the data in a particular direction.

A more rigorous method, as shown by Turcotte [1992], may be used to determine the two dimensional spectra of a surface. Consider a NxN data matrix of length L representing equally spaced measured data values from a two dimensional surface. The N² data points are denoted by h_nm with (n,m) indicating the location of the data points in the x and y directions respectively. The following steps, as outlined by Turcotte (1992), are carried out to determine the two dimensional spectra:

Step 1: Compute the two dimensional DFT and deternine the complex coefficients, H_st as follows

where s=0,1,...N-1 and t=0,1,...N-1 are the transforms in the x and y directions respectively.

Step2: Assign each coefficient H_st an equivalent radial number, r, using

Step 3: Compute the 2D power spectral density, S_2j, for each radial wave number, k_j, using

where N_j is the number of coefficients that satisfy j<r<j+1 and the summation is carried out over all H_st in this range. Finally, the b parameter is estimated from a linear fit to the slope of a log-log plot of power spectral density versus radial wave number. The relationship between power spectral density and wavenumber for scale invariant data is

We use the method outlined above to estimate the two dimensional power spectra of our seismic data surfaces extracted from within the fluvial and deltaic systems.

Estimating the parameter b:
The parameter b is determined from a least squares linear fit to the slope of a log-log plot of power spectral density, S(k) or P_CAP^(k), versus wavenumber, k. The sampling frequency, 1/delta, and the length of the data record, N, determine the useful portion of the log-log plot to be used in determining the parameter b. Thus, the linear fit is determined over a range of frequencies (or wave numbers in the case of spatial data) for which the spectral estimate is highly accurate. Aliasing problems tend to occur at the high wavenumber end of the spectrum and a first order estimate for determining the high wavenumber cut-off is to discard wavenumbers above one-half the Nyquist frequency (1/2delta). The high wavenumber cut-off is thus 1/4delta. The low wavenumber cut-off depends on the length, N, of the data record and the longest lag for which the autocorrelation is estimated. We use low wavenumber cut-off values of 4/Ndelta and 1/2pdelta for the modified periodogram and Capon estimator methods respectively [Austin et al., 1994].