A generation of geophysicists in the petroleum industry has successfully utilized seismic attributes for interpretation of the porous zones related to hydrocarbon accumulations. While interpretting permeability using these techniquesis somewhat more difficult, on a local field scale, some success has been realized. These same techniques could provide valuable input to hydrogeologic studies by contributing an estimate of the distribution parameters of correlation length and fractal dimension of hydraulic conductivity.

This study provides a relationship between seismic reflectivity and hydraulic conductivityfrom the literature. A non-traditional seismic seismic survey has been acquired in a manner allowing for future validation of the interpretation. This survey is used as input to geostatistical methods providing the correlation length and fractal dimension for the formation under study. These results are presented along with a critique of the methods utilized for these interpretations.

Prediction of the stochastic hydrogeologic properties of correlation length and fractal dimension of an aquifer's hydraulic conductivity are valuable for hydrogeologic studies that utilize stochastic modeling techniques. While hydraulic conductivity can be calculated from permeability that has been measured directly from core samples obtained from wells or calculated from a slug test, interpretation of the distribution of the hydraulic conductivity and the associated stochasic properties of the formation(s) of interest is more problematic. The utilization of a remote sensing tool that responds to a physical property related to hydraulic conductivity should supply insight to these properties. The modern high resolution seismic survey is an inexpensive method that can supply the data for this purpose.

The relationship between seismic reflectivity and hydraulic conductivity is established and methods of extracting correlation length and fractal dimension from a seismic survey are presented. The seismic survey acquired for this work has been designed with a nontraditional geometry in a manner such that the interpreted stochastic parameters may be confirmed with future permeameter studies.

Seismic data has long been used in the field development stage of the petroleum industry. In practice, seismic interpretation has successfully contributed to the identification and exploitation of producing zones in oil and gas fields because seismic methods respond to a property (seismic velocity) that within these producing zones is dependent on porosity. Great confidence in this method has been gained in the petroleum industry from the number of successful applications over last several decades of use.

While target zones in the petroleum industry are typically at considerably greater depths than zones of interest for hyrogeologic studies, the parameters of interest are similar. Both pursuits benefit from knowledge of permeability and porosity values as well as distribution of these properties. While several generations of geophysicists have applied these methods in the petroleum industry, these methods have yet to be applied to hydrogeologic studies. Utilization of shallow, high resolution seismic methods for structural purposes has been incorporated in environmental and engineering geophysics with good success since the early 1980's. Steeples et al. (1986) and Knapp and Steeples (1986a and b) document the progression of such a project from acquisition through interpretation of standard common midpoint (CMP) data for such a survey.Interpretation of seismic attributes such as that done in the petroleum industry will advance the use of the seismic method for surveys such as these. The objective of this study is a preliminary step towards this goal. The goal is really three-fold: 1) establish from the literature a link between seismic response and hydraulic conductivity, 2) design a seismic survey at a site where future evaluation of the results with permeameter studies is possible and 3) acquire a seismic survey, process and interpret the seismic reflectivity data for correlation length and fractal dimension. Accomplishing this will result in an interpretation of stochastic parameters useful in hydrogeologic modeling studies. Future permeameter studies will be necessary to evaluate these predictions.

Acquisition of seismic data in a method such that a future permeameter survey can fully evaluate the results required a ontraditional geometry of the area to be surveyed. Ideal geometry included a relatively flat surface for source and receiver locations and a nearly straight face for a distance longer than the correlation length of the formation. Additionally, good access to the site was desirable as transportation of equipment over long distances would hamper acquisition. The Northwestern Nevada and Northeastern California area has many sites satisfying some of these attributes to varying degrees. Over a period of five months in 1995, a number of these sites were evaluated for suitability for this study. An abandoned open pit diatomite mine possessing a high degree of these desired features was chosen for this study (see Fig. 1 and Fig. 2). The site chosen, near Hazen, Nev. also had a relatively short commute distance of some 40 miles (see Fig. 3).

Design of the seismic survey had to conform to the geometry of the mine bench where data was acquired. The bench chosen could accommodate a survey to include one line parallel to the face of the mine and three lines transverse to the face and parallel line. Two days were required for the design, location and elevation survey for acquisition.

Test records were acquired in the field on the first day of acquisition to determine recording instrument settings, optimum source parameters and minimum offset distance. After determining acceptable field acquisition parameters, three full days of acquisition resulted in the raw field data for four high resolution seismic lines that are the focus of this study. Data acquisition work was limited to early morning through early afternoon hours due to the high afternoon temperatures typical of the desert climate where the survey was acquired. Downloading the data from the recording instruments and the laptop computer used for temporary storage was accomplished during the afternoon periods.

Raw field data was processed into brute stacks, two dimensional and three dimensional migrated sections over a seven month period. Power spectra and semivariograms were calculated for all four lines and a three dimensional migrated section that was constructed at a face location. This processing sequence spanned two months. The output from the power spectra and semivariogram processing has been interpreted for correlation length and fractal dimension. That interpretation is presented here.

Hydrogeologic studies are often the major component of programs required for the remediation of contaminated groundwater. Location of the contaminant(s) of concern within the aquifer of interest, movement of groundwater and contaminants and the distribution of the concentration of the contaminant with time is of primary interest of these programs. Greater understanding of the hydrogeologic parameters of the aquifer of interest will lead to better design of these remediation procedures.

The determination of a contaminant's concentration within an aquifer throughout time and space is predicted by the stochastic transport governing equation as expressed by Gelhar et al. (1979):

where A is the stochastic macrodispersivity in equation 20 of Gelhar et al., alpha_L is local longitudinal dispersivity (in L units), U is the mean velocity (in L/T) units and xi is a moving coordinate system. Using this governing equation, the distribution of concentrations of a contaminant throughout time for a two dimensional cross section in a study area can be determined.

Recent advances of computer speed and memory have permitted successful implementation of stochastic modeling programs for analysis of contaminant problems. With these advances, field scale stratified aquifer models can be constructed. To construct these models, stochastic macrodispersivity values for these models must be determined.Gelhar et al. (1979) provides us with a stochastic method to determine dispersivity in two dimensional cross sections of stratified aquifers. At distances greater than the correlation length of the formation of interest, Gelhar develops a value for longitudinal dispersivity determined by use of:

where A is asymptotic macrodispersivity, sigma is the standard deviation of the log of hydraulic conductivity (in L/T units), l is correlation length, K bar is the average of the formation's hydraulic conductivity and alpha_T is the transverse dispersivity.Upon determining this value, it can be used in equation 1.

Determination of asymptotic macro-dispersivity from equation 2 necessitates not only knowledge of absolute values of hydraulic conductivity of the formation of study, but also knowledge of its correlation length for hydraulic conductivity. While hydraulic conductivity itself can be measured from cores, outcrops or other methods utilizing direct physical contact with the formation, distribution characteristics such as correlation length must rely on remote sensing methods. While calculation of the average and standard deviation of hydraulic conductivity relies on knowledge of absolute values of this parameter, distribution parameters need only utilize relative values. It is this characteristic of correlation length and fractal dimension that allows extraction of these parameters directly from seismic data throughout the entire volume under study.

In addition to correlation length the statistical distribution parameter of fractal dimension can be extracted from seismic data. Some models for fractal theories of transport in aquifers utilize fractal dimension for macro-dispersivity given volume from an aquifer. Many geologic properties have been shown to be accurately modeled by fractional Brownian motion. Such distributions are characterized by fractal dimension. These methods may beneficially effect extraction/mitigation techniques of remediating contaminated groundwater by estimating the time needed for extraction of the contaminated volume of water. Recent research such as Wheatcraft and Tyler (1988) and Zhan and Wheatcraft (1996) demonstrate that fractal dimension is a measure of tortuosity of the aquifer. Knowledge of this parameter will lead to better evaluations for the solution to the contaminant transport problems.

Seismic waves respond to the physical properties of the volume of material surveyed in three basic wave attributes: amplitude, frequency and phase. The amplitude of a seismic wave is a response to the reflectivity of the formation boundaries under investigation and is the attribute examined for this study. Reflectivity has long been used to extract seismic acoustic impedances and velocities from seismic reflections for depth conversion, formation thicknesses and other applications. Nearly a generation of petroleum geophysicists have successfully utilized seismic acoustic impedances and velocities extracted from reflection coefficients to predict porosity zones in oil and gas reservoirs at considerably greater depths than the typical hydrogeologic study. The analysis of high resolution, shallow seismic data in a similar fashion will produce similar results.

Recent research has established a good correspondence between correlation length of seismic reflectivities extracted from standard CMP processing of deep crustal reflectors and associated correlation length of seismic velocities collected from outcrops of the formation under study (Levander et al., 1994). Pullammanappallil et al. (1996) advances this thread of research by extracting correlation length and fractal dimension from the power spectra of deep crustal reflectors. These research projects establish the use of standard CMP seismic data for determination of seismic velocity stochastic parameters.

As early as 1950, petroleum industry researchers recognized that within lithologically consistent stratigraphic units (all clastics or all carbonates) porosity could be used to predict permeability within field scale areas of interest (Wylie and Rose, 1950). These methods have been applied with some success in the petroleum industry (Schlumberger, 1972a,b). This relationship between porosity and permeability leads to the derivation of hydraulic conductivity and supports the use of correlation length of seismic reflectivity for prediction of correlation length of hydraulic conductivity. The process proceeds from seismic reflectivities to estimating acoustic impedances, which estimate seismic velocities, which estimate porosities. On a field scale study within a single formation, porosities can predict permeablities which can be converted to hydraulic conductivities. By establishing this relationship between seismic reflectivity and hydraulic conductivity, the relationship between correlation lengths of these parameters is also established.

Dobrin (1960) defines the seismic reflection coefficient (reflectivity) as:

where rho₁ and rho₂ are the densities of the formations (in M/L³ units), v₁ and v₂ are the respective seismic velocities (in L/T units) and r is reflectivity. Lindseth (1979) derives an empirical relationship for density as a function of velocity where:

where density is measured in gm/cm³ and velocity is in ft/sec. (The incompatible units are kept here to keep equation 4 in the same form as the original paper.) Lindseth combines equations 3 and 4 producing an inversion section displaying seismic velocities.

Velocity data derived in this manner has been successfully utilized for porosity analysis in the petroleum industry for two decades. A common method for determining porous zones from sonic logs is included in Schlumberger (1972a) as:

where phi is porosity, v_o is observed seismic velocity, v_m is matrix velocity and v_f is velocity of the pore fluid (velocities are all in L/T units). Porosity extracted in this manner can be used on a field scale study for a single formation as a predictor of permeability (Wylie and Rose, 1950 and Schlumberger (1972a,b). Schlumberger presents the relationship between permeability and porosity as:

with k being permeability (in L² units), phi porosity, s_w irreducible water saturation and C is an empirical constant after Schlumberger. Permeability can then be converted into hydraulic conductivity using:

where rho is density, mu is viscosity (in M/LT units), g is gravitational acceleration (in L/T² units), k is permeability and K is hydraulic conductivity. Combining equations 5,6 and 7 with the seismic velocity derived from equations 3 and 4 hydraulic conductivity can be calculated with:

where M is a constant determined by:

This development establishes that within an aquifer an absolute value for hydraulic conductivity can be determined dependant only on parameters that can be expected to be reasonably consistent throughout the aquifer of study (viscosity and density of the fluid, irreducible water saturations, formation constant, fluid velocity and matrix velocity) and one parameter calculated from eqns. 3 and 4 (seismic velocity).

While direct prediction of absolute values of hydraulic conductivity from seismic data may be difficult due to the presence of the empirical constant of equations 5 and 8, it should be noted that values for M and v_m can be expected to be reasonably constant for any given aquifer on a field sized scale. The absence of great horizontal variation of these properties within the aquifer of interest allows for extraction of meaningful values for both correlation length and fractal dimension of an aquifer from reflection seismic data. These interpreted stochastic parameters are related to the corresponding stochastic parameters of the hydraulic conductivity of the formation and can be extracted by use of semivariograms. The value of correlation length thus interpreted for can then be used in determining dispersivity as in equation 2. Fractal dimension interpreted in this manner can be used in determining tortuosity and pump time for remediation methods.

Statistical analysis of high resolution seismic sections used for interpreting correlation length and fractal dimension of reflectors can be initiated using two processes, the power spectrum and the semivariogram. These processes, while both contributing insight to an interpretation of these distribution properties, differ in method of calculation and application. These differences can be used to enhance the interpretation of these stochastic parameters.

Power spectra extracted by utilizing the fast fourier transform (FFT) method have been a common component of seismic processing used for analysis of frequencies since the mid 1960's. The process transforms seismic data from the time domain to the frequency domain. Press et al(1986) expresses the power at a given frequency as:

where k represents the k^th component of the corresponding transform. The result can be plotted as a log of total power vs. log of frequency. This same method can be adapted to the horizontal direction of a seismic section resulting in log of power plotted against log of spatial frequency. The resulting power spectrum will contain seismic responses to the input signal and lithology. These power spectra can then be interpreted for the correlation length and fractal dimension of the seismic response to the lithology for both the vertical and horizontal components. The peaks observed on the power spectra correspond to the wavelengths of highest correlation for both seismic signal and lithologic response. The lower responses of these plots correspond to spatial frequencies with little or no correlation. Distances beyond these peaks correspond to distances beyond which white noise is prevalent. For these distances stochastic models are valid and the correlation length interpreted can be used in equation 2.

A semivariogram is a plot of the sum of variance of points in a given field and its neighboring points against a distance between the point under consideration and the points within a specified distance window of this point (called lag distance). The semivariogram can be calculated as in Carr (1995) by:

where Z(x_i) is the value of the parameter under consideration and Z(x_i+h) is the value of the parameter at the neighboring point(s) within the lag distance window. An alternative estimator of the semivariogram using a moving window can also be used. The plot of the semivariogram estimator after Li and Lake (1994) results in a what is described as an estimator that is "unbiased, robust and resistant to contamination". This estimator is more stable and the results are more easily interpreted than the classical semivariogram of equation 11, the moving window calculation acts as a filter for the plot. This estimator is calculated by:

where n is the number of data values in the entire field, D_i.h is the index set of data values in the moving window (of size h centered at block i), excluding point x_i , and m is the number of data values in D_i,h.

A linear plot of either the classical semivariogram or the moving window estimator against the lag distance for most geologic attributes results in a graph with characteristics of a positive slope increasing to a sill characterized ideally by a horizontal line. The value of this sill can be interpreted as the variance of the dataset. The lag distance where this value is reached, called the range, can be interpreted as the correlation length used in equation 2. Beyond this distance there is no further correlation between the values of the points and stochastic models are valid for the parameter under study.

Correlation length and fractal dimension can be interpreted from both power spectra and semivariograms. Carr and Myers (1984) and Carr (1997) characterizes most semivariogram models of geological data to be of the spherical model (see Carr, 1995). This model reaches a well defined sill at a well defined range. The range of the semivariogram indicates the distance at which greatest variance is reached for the parameter under study. The correlation length interpreted from this model is the range of the semivariogram. Power spectra emphasize correlation of data as opposed to the emphasis semivariograms place on variance. Low values of correlation interpreted from power spectra in the wavelengths related to lithologic variations can be interpreted as correlation length of the formation when no further peaks are indicated at lower spatial frequencies corresponding to longer distances. This distance may also be interpreted as the correlation length of the data field. Utilization and comparisons of both of these methods serve to enhance the interpretation of this parameter.

Fractal dimension of a field of data can be interpreted from the slope of the power spectrum. Carr (1995) uses:

where b is the absolute value of the slope of the power spectrum. For this purpose, this should be applied only to portions of the spectrum that correspond to spatial frequencies relating to lithologic variation. The slope of the semivariogram can also be used to determine fractal dimension by:

after Carr (1995) where S is the slope of the log of the semivariogram plotted against log of the lag distance. Here again, care must be given to apply this only to portions of the semivariogram corresponding to lithologic variations. Utilizing both methods here raises confidence in the interpretation when the results are similar.

In July, 1996, four seismic lines were acquired on a bench of Eagle Picher's Section 8 diatomite mine located in sec. 8 of T 19N R 26E some 3 miles southwest of Hazen, Nevada in Churchill County and about 40 miles east of Reno. The diatomite is deposited in a Tertiary sedimentary sequence described by Moore and Archibold (1969) as "Lacustrine and fluvatile sediments. Sandstone, shale, marl, diatomite, limestone and calcareous tufa. ... Includes Truckee Formation, Aldrich Station, Coal Valley, and Morgan Ranch Formations of Axelrod (1956)." Figures 4 and 5 are closeup pictures of the face and top of the bench respectively. The bedding varies on a scale of millimeters to tens of centimeters with thin interbedded darker colored sands and shales.

A Bison 9000 12 channel recording system on loan from the Desert Research Institute was used to record nominally 24 fold data. The source was a 3 lb. single jack hammer. Ten strikes of the hammer against a hardened steel milling ball were summed for each record. Twelve foreshot and twelve backshot records at offsets varying from 36 to 6 feet were acquired for each geophone setup to achieve 24 fold coverage. Recording filters were set as a 64 Hz low cut filter and out for the high cut. Records were recorded with a .25 millisecond sampling rate resulting in a 2000 Hz Nyquist frequency. Recording time was 500 milliseconds.

Four lines were acquired, one for a length of 144 feet parallel to the exposed face of the pit at a distance of approximately 50 feet from the face and three lines transverse to the face and first line, tieing the first line and ending at the face (see Fig. 6). The line parallel to the face utilized a 2 foot geophone spacing with a 2 foot source spacing. Geophone setups were recorded with twelve foreshots and twelve backshots then "moved along" six stations where the process was repeated. The resulting 24 fold seismic line has a CMP spacing of 1 foot (0.3048 m). Transverse lines were acquired using 1 foot geophone spacing and a 2 foot source spacing. The result is 24 fold coverage with a CMP spacing of 1/2 foot (0.1524 m).

Acquisition required three days of work. Weather was hot, dry and sunny with no wind. Very little noise was observed on the raw data records and in processing one record and a single trace for each record of one setup were all that were discarded as unusable for a total of 28 of the 5460 traces recorded.

Standard CMP processing as described by Dobrin and Savit (1988, pg. 259) was applied to all four lines. Brute stacks at a series of 19 stacking velocities ranging from 350 m/sec to 11,000 m/sec were generated and used to determine the reflection from the face as the higher velocities enhanced the face reflector whereas normal reflections from bedding deteriorated at the high velocities. With nearly horizontal bedding (see Figs. 4 and 5), the only event demonstrating apparent dip in the brute stacks of the transverse lines (Figs. 7, 8 and 9) is the face reflector. Both two dimensional and three dimensional migrations were generated from the data set. The three dimensional migration used was the three dimensional Kirchoff summation method as described in Claerbout (1992 , pg. 106).Figures 10 through 14 (fig.10, fig.11, fig.12, fig.13, fig.14) are migrated sections placed at a position near the face of the pit and lines of acquisition respectively. These lines represent only the upper 150 feet of data as reflections from the face are entirely positioned by that depth.

Standard CMP stacking of steeply dipping events can result in attenuation or even destructive interference of the signal of interest (Gibson and Levander, 1988). This dip filtering effect must be avoided for power spectra and semivariograms to be extracted from any cross sectional image. Avoiding this dip filtering effect can be accomplished by use of pre stack migration techniques such as the three dimensional Kirchoff migration method. This technique properly images points within thedata volume and is applied to both the horizontal and vertical components in the same manner resulting in no distortion of the imaging or corresponding statistical parameters extracted here. This migration technique does, however, introduce some artifacts which affect the edges of the migrated sections and produces spurious arcs due to badly gained traces and noise. These artifacts can readily be identified and are not a hinderance to the extraction of the statistical parameters of interest.

Louie et al. (1988) describe the three dimensional Kirchoff migration method utilized here as particularly useful for imaging vertical or near vertical events such as the face of the pit studied here. The face line was properly positioned using this method and many features of the face have been imaged with the migration. Comparing a picture of the face (fig. 4) to the three dimensional migrated face line shows the resolving power of three dimensional migration. Structural features of the face such as the fracture zones, ledges and some bedding has been imaged. Migration of the transverse lines images the face in the proper position for these lines confirming the value of this method. Proper placement of geologic features and lithologic variations utilizing this method lead to a more confident interpretation of the distribution parameters of interest for this study.

Semivariograms were calculated for all migrated sections. A classical semivariogram and the moving window semivariance estimator after Li and Lake as presented previously were used for calculating semivariograms. Semivariograms were calculated in both the horizontal and vertical directions for all of the migrated sections. Figures 15 through 19 (fig.15, fig.16, fig.17, fig.18, fig.19) are generated for the horizontal component. Figures 20 through 24 (fig.20, fig.21, fig.22, fig.23, fig.24) are the semivariograms for the vertical component. The classical and Li and Lake estimator are plotted together on all plots. The filtering characteristic of the Li and Lake estimator proved to be more interpretable for this dataset. These plots were replotted on log-log scales (Figs. 25 through 29) (fig.25, fig.26, fig.27, fig.28, fig.29) to determine slope for fractal dimension interpretation. In addition, values for the shorter lag distances corresponding to seismic signal variance were deleted. This technique produces plots with slopes more easily interpreted than the classical method.

Power spectra were extracted from all depth migrated sections using the fast fourier transform method as above. Spectra were generated in both the horizontal and vertical directions for these sections (Figs. 30 through 34) (fig.30, fig.31, fig.32, fig.33, fig.34). The output contains peak values at the spatial frequencies corresponding to the highest correlation values of the data as they are correlated with known frequencies. A log-log plot of this correlation against the spatial frequency used for the correlation is then generated. Here the log of the spatial frequency has been converted to their corresponding wavelengths for ease of interpretation. Note that the highest frequencies plotted on the power spectra (corresponding to the shortest wavelengths) are controlled by the Nyquist frequency of the migrated sections (2 foot or 0.6096 m). High power values represent high correlation and low values, low correlation. The value where power decreases in the portion of the power spectrum corresponding to lithologic changes can be interpreted as the correlation length of the data field. While we can interpret a correlation length from line A and the face line, the power spectra generated for the shorter lines demonstrate a continuous increase with wavelength throughout the portion responding to lithology. A reliable correlation length is not interpretable from these lines due to their length and lack of datapoints generated with this method.

The slope of these power spectra can be used to determine fractal dimension. Again the low frequency portion of the power spectra that corresponds to the lithologic data is sparsely controlled for the shorter lines leading to interpretations of little value utilizing solely this method. However, these lines where few datapoints control the slopes used for the interpretation can be used to confirm values obtained utilizing other methods.

Power spectra and semivariograms (both classical and Li and Lake's estimator) have been generated from the migrated data of the four lines of acquisition and the face line. An interpretation of the semivariograms (in both linear and log-log plots) and power spectra are presented here. The interpretation of the power spectra is included here merely for completeness. The lack of control in the bandwidth corresponding to lithologic response displayed by the power spectra contributes to difficult interpretations from this method. This lack of control is discussed in detail in the comparison of methods section.

Figures 15 through 19 (fig.15, fig.16, fig.17, fig.18, fig.19)are plots of both the classical and moving window estimators of the semivariogram. The shorter plots (Figs. 17 through 19) were extended to the length of the plots of the longer lines by adding one point on the moving window estimator with a gamma value of 1.0 at the longest lag distance of the longer lines. This results in all semivariograms plotting at the same horizontal scale. Figures 20 through 24 show the vertical component and are included here for completeness. These figures show little more than the seismic signal contribution. The single steep slope seen on the plot of the vertical component is characteristic of the seismic signal with a peak wavelength of 1.4 m. This corresponds to point A plotted on the face power spectra (Fig. 30) which is the high frequency contribution of the seismic signal. The horizontal component semivariograms contribute more to the interpretation. These semivariograms show good development of two distinct slopes before reaching a sill. This is indicative of variance contributions from two sources, the seismic signal and the lithologic variation of the formation. The first steep slope and change in slope corresponds well to the vertical component. This is variation due to the seismic signal. The second, gentler slope from a distance of approximately 3.1 m to the range is the contribution of lithologic variation to the semivariogram. This second slope is useful for interpretation of correlation length and fractal dimension for the lithology sampled.

The semivariograms extracted from the data analyzed here exhibit characteristics of the spherical model. Interpreting the semivariograms with this model allows for a range (and corresponding correlation length) relating to a high percent of variation near the fully developed sill. The characteristic shape of these plots displays an asymptote to the sill as lag distance increases to 100 percent variance. As the plot approaches the sill asymptotically, a change in slope is apparent on all the horizontal semivariograms. Interpreting correlation length at the lag distance corresponding to 100 percent of the variance results in lengths of perhaps longer than desirable. Picking the change in slope returns a result that is perhaps shorter than desirable. The lag distance interpreted as the range for this study is a point halfway between this change in slope (point A in Figures 15 through 18) and the distance corresponding to the point where the sill of the semivariogram becomes well developed (point B in Figures 15 through 18). This interpreted correlation length is noted with an asterisk in these figures. The horizontal semivariogram extracted from line D (Figure 19) does not demonstrate a long asymptotic approach to the sill. A correlation length can therefore be directly interpreted as the lag distance corresponding to the point of the well developed sill (noted by the asterisk). This interpretation results in correlation lengths of approximately 18 m, 18 m, 13 m, 10 m and 15 m respectively for the face line, and lines A through D. The north-south lines (face line and line A) demonstrate good agreement with a correlation length of 18 m from the spherical model. This can be interpreted as the northing component for this formation. The east-west lines (lines B, C and D) also show a good agreement with interpreted correlation length between 10 m and 15 m. This range can be interpreted as the easting component of the formation.

An alternative interpretation utilizing the method of Gelhar (1993, pg. 311) and an exponential model for the semivariogram utilizes a lag distance corresponding to a gamma value of:

(where sigma² is the sill of the semivariogram) is determined from the semivariogram and used as the correlation length of the sampled stratigraphic sequence. This corresponds to 63.2% of the variation of the data sample. The variation due only to the lithologic contribution of the semivariogram is utilized for this analysis. Therefore the lag distance is interpreted from the second slope of the semivariogram. The correlation lengths interpreted using this method are approximately 11 m, 11 m, 9 m, 7 m, and 10 m for the face line and lines A through D respectively.

A method of extracting fractal dimension from Carr (1995) utilizes the slope of a log-log plot of the semivariogram and the relationship:

where S is the slope of the semivariogram from this plot. Figures 25 through 29 (fig.25, fig.26, fig.27, fig.28, fig.29) are log -log plots of the face line and lines A through D. Data from the first 3.1 m have been deleted for ease of interpretation of the lithologic contribution to these plots. For this interpretation, a line has been visually placed on the linear portion of the plots (i.e. over a log of lag distance range of from 0.5 to 1.2 on Figure 25) . The slopes were then measured and fractal dimension determined from equation 16. The slopes of these plots are 0.12, 0.12, 0.26, 0.27 and 0.30 respectively. The resulting fractal dimensions are 1.94 for both north-south lines and 1.87, 1.86 and 1.85 for the east-west lines.

The vertical and horizontal power spectra are plotted together in Figures 30 through 34 (fig.30, fig.31, fig.32, fig.33, fig.34) so the contribution of the different components can readily be identified. Note that the frequency band of the vertical component is entirely in the range of the frequencies of the seismic wavelet. The low cut filter used in recording was 64 Hz (relating to a wavelength of approximately 5.5 m). The data was recorded with a 0.25 ms. sampling rate for a Nyquist frequency of 2000 Hz. The peak frequency (175 Hz.) is well within that range.

Both the vertical and horizontal components of the power spectra contain information from both the seismic wavelet and lithology. However, even for the vertical power spectrum exhibiting the broadest bandwidth (the Face line seen in Figure 30), it can be seen that the bandwidth is limited to the 2 m to 8 m range which corresponds to the frequencies of the seismic wavelet. The vertical response to the seismic wavelet overpowers what little seismic response there is to the thin bedded lithology with little or no vertical velocity contrasts in this area. A similar response in this bandwidth is illustrated by the horizontal components of the power spectra. Here again the response to the seismic wavelet is overpowering lithologic responses at the shorter wavelengths. The horizontal component does display a good second peak corresponding to the longer wavelength variation in lithology. It is this variation that can be utilized to extract the correlation lengths used in determining macrodispersivity for this formation.

At the longer wavelengths (corresponding to the lower frequencies), the horizontal and vertical components of the power spectra separate. The horizontal component demonstrates a second maximum. This maximum develops completely for the longer lines (the face line, line A and line D) falling away from peak power. This indicates that a correlation length is reached for these lines. This peak is reached at lengths greater than 8m and is attributable to lithologic properties. It is this range of wavelengths that can be utilized for interpretation for the stochastic properties of interest here.

Correlation length can be interpreted from power spectra by determining the length where correlation falls off from the peak corresponding to the lithologic properties. This peak can be seen on power spectra generated from the longer lines at approximately 20 m. The correlation length can therefore be interpreted as between 20 and 40 meters as the next control point of control on these plots is 40 m. This value, while controlled by fewer datapoints than ideal, does confirm the interpreted value obtained from semivariograms. Note the sparsity of data used for this interpretation from the power spectra. This lack of data is due to fast fourier transform methods of calculation. Further discussion of this drawback is included when the methods of interpretation are compared.

Fractal dimension can be interpreted from the slope of the power spectra in the portion of the plot corresponding to lithologic variation. However, interpretation of fractal dimension independently from the power spectra extracted from this dataset is difficult due to the lack of control in the bandwidth corresponding to lithologic variation. Lines with slopes corresponding to the fractal dimension extracted from the semivariograms have been plotted on the power spectra through the portion corresponding to lithologic variation. These slopes of these lines demonstrate that the fractal dimensions from the semivariograms are similar to values that could be interpreted from the power spectra. These plots only serve to confirm the interpretations from the semivariograms over limited bands of the power spectra.

The good agreement between the stochastic parameters as observed with this dataset raises the confidence placed in the values obtained by this interpretation. While similar values for both correlation length and fractal dimension are obtained for each component minor differences are observed. Good agreement between the results of the interpretation can be seen in the northing component of these parameters. While the east- west lines also demonstrate high levels of agreement of these parameters between lines, the interpreted correlation length of between 10 m and 15 m being nearly the length of the shortest line (line C at 15.2 m) is a concern for the values extracted for the easting component from this line. The other two lines however demonstrate good sill development and are in good agreement on the stochastic parameters interpreted here.

While the power spectra and semivariograms presented here have been generated from the amplitudes of migrated seismic sections, previous studies (Levander et al, 1994 ; Pullammanappallil et al., 1996) have demonstrated a relationship between these stochastic parameters extracted from seismic reflectivity and corresponding velocity data. Similar correspondence between these same stochastic parameters extracted from velocity fields and the calculated hydraulic conductivity field should also be established. A random seismic velocity dataset was evaluated to study this relationship. The following section is a discussion of this study.

Equation 8 presents hydraulic conductivity as a function of seismic velocity. This equation was used to convert a dataset of 300 random numbers representing velocity values of a clean sandstone with porosities of between approximately 2% and 30%. A value of M was chosen so that these velocities would be converted to realistic hydraulic conductivity values from Bear (1979, pg. 68). These values were structured as if they were a 5 foot by 60 foot (1.5 m by 18.3 m) volume of a formation. Semivariograms were calculated for both the input velocities and the resulting hydraulic conductivities. Figures 35 and 36 are the plots of the semivariograms constructed from Li and Lake's moving window estimator. The semivariograms of the random sample demonstrate good similarities along the entire plot as would be expected from totally random samples. A second dataset was constructed from these 300 random velocities with the first sixty samples unchanged, a constant value of 500 ft/sec was subtracted from the second sixty values, 1000 ft/sec was subtracted from the third sixty samples, 1500 ft/sec from the fourth, and 2000 ft/sec from the fifth sample. The result is a major slowdown of velocity every 12 feet (or 3.66 m) producing a correlation length in the dataset. Again semivariograms were constructed for both the velocity values and the hydraulic conductivities. Figures 37 and 38 are the plots of these results. There is good agreement of these semivariograms for the range. Both plots exhibit a range value of approximately 1.5 m. While many features are similar, such as the range of the semivariogram there are some differences. The major dip in values of gamma at lag distances just longer than the range is considerably more significant for the calculated values of hydraulic conductivity than the velocity data. Log-log plots of the semivariograms were produced for facilitating fractal dimension interpretation. These plots (Figures 39 and 40) somewhat subdue the difference in the minimum seen on the linear plots. Slopes of these plots for the portion of the semivariogram where log of the lag distances are shorter than the range are .098 and .184 respectively resulting in fractal dimensions of 1.95 and 1.91 for the velocity data and hydraulic conductivities. While this difference does not appear large in an absolute sense, the relative difference in slope is considerable as the slope of the semivariogram for the hydraulic conductivity is nearly twice the slope of the semivariogram from the velocity data. It is not readily apparent from equation 8 why the slopes would differ from each other by a factor of two.

It should also be noted here that the slope under consideration for this parameter is controlled by only five of sixty points. Additional studies may even determine that this difference is statistically insignificant. Future investigations into the relationships between the slopes of semivariograms of hydraulic conductivity data calculated from velocity data will be necessary for full utilization of fractal dimension extracted in this fashion.

Figure 41 is a picture of the face in a position nearly duplicating the migrated face section. The vertical boards are place for scale and location. They are 8 feet in height with alternating blue and orange 1 foot sections and are placed at the tie points of the transverse lines. This picture was taken with a digital camera with a 640 by 480 pixel resolution. The picture was then cropped in photoshop to delete all data except the face. The cropped image was then resampled to nearly match the migrated seismic section sampling. The result is a 107 foot by 24 foot (32.6 m by 11.0 m) visual sample of the face. Power spectra and semivariograms were then calculated for this image. Figure 42 is the power spectra, figures 43 and 44 are the horizontal and vertical components of the semivariogram respectively and figure 45 is the log-log plot of the semivariogram of the horizontal component. Applying the techniques described above results in a 17 m correlation length interpreted from the semivariogram and fractal dimension of 1.77. The slope corresponding to this fractal dimension has been plotted on the power spectra extracted from the photo in a similar manner as has been done with the seismic data. This plot again shows that the fractal dimension interpreted from the semivariogram can be confirmed from the power spectra.

A summary of the correlation lengths and fractal dimensions interpreted from all of the lines and the photo are presented in table 1 for ease of comparison. This interpretation is from the semivariogram only as this method is more reliable as previously discussed. From this comparison it can be seen that the correlation lengths extracted from the north-south datasets (the photo, Face Line and Line A) demonstrate excellent agreement (within 1 meter). Correlation lengths extracted from the two longer east-west datasets (Line B and Line D) also show an excellent agreement (within 2 meters). Correlation length from the shortest line (line C) may be affected by its length as previously noted but is still within good agreement for that component. Fractal dimensions extracted from the north-south seismic lines also demonstrate excellent agreement. Considerable difference between the seismic data and the visual data is noted from this comparison. While 0.17 is considered a significant difference in fractal dimension (see Wheatcraft and Tyler (1988) for a demonstration if a fractal dimension difference of 0.2 and its effect on tortuosity) it can be expected that seismic data will inherently display higher fractal dimension than visual data due to noise (see discussion below). The comparison of the fractal dimensions extracted from the east-west seismic data shows excellent agreement with a difference of only 0.02 from the biggest to the smallest value extracted.

Horizontal ............................Horizontal

Correlation Length................ Fractal Dimension

Photo ....................17 m .....................................1.77

Face Line ..............18 m .....................................1.94

Line A ...................18 m .....................................1.94

Line B ....................13 m.................................... 1.87

Line C................... 10 m .....................................1.86

Line D................... 15 m .....................................1.85

The excellent agreement between the correlation length extracted from the photo and the seismic data demonstrates the value of using seismic data for interpreting this property and utilizing it for other lithologic parameters. While the fractal dimensions interpreted in a similar manner result in higher values from the seismic data than the visual spectrum, the introduction of white noise to the compared datasets should be considered. While the data from the visual spectrum has little contamination of white noise, the seismic data has a considerable level of noise for even the best of acquisition techniques. This noise will result in an overall increase in fractal dimension for the seismic data. The resulting interpretation of the fractal dimensions from these data should therefore be considered as an upper limit for this property. It should also be noted that the noise level inherent in the seismic data will be the same for all three components (easting, northing and vertical) will be unaffected. These relative differences between interpreted fractal dimensions for the different components can be useful in interpreting macrodispersivity.

Correlation length of stratigraphic zones can be interpreted by use of semivariograms and power spectra using the methods presented here. The ranges of the semivariograms are easily interpreted from the plots acquired from this dataset resulting in correlation length for the component under study. Here there is good agreement between the correlation lengths acquired from the lines in the same direction and the interpretation benefits from such agreement. Interpretation of correlation length from power spectra is more problematic. The logarithmic result of the power spectra gives little control over the portion responding to the lithologic variation of this dataset. The longest lines (line A and the face line) return only 9 datapoints that record information at lengths that can be attributed solely to lithology. This compares to 131 points controlling the lithologic portion of the semivariogram. While correlation length is difficult to interpret from power spectra, a result similar to that interpreted from the semivariograms can be interpreted. This agreement from an independent method serves to validate the interpretation of correlation length. Power spectra results merely confirm the results of the semivariogram interpretation and demonstrate the ability to differentiate between seismic signal and lithologic variation in the dataset.

Fractal dimension can also be interpreted independently from semivariograms and power spectra. The useful distance of analysis for this dataset falls in the range of 5.5 m to 45 m (the wavelength corresponding to the low cutoff frequency used in recording the data to the full length of the longest line). This useful length is a function of the equipment and field procedure used. The semivariogram procedure returns a linear plot over the length of the line whereas the power spectrum is a logarithmic result. Again this results in the lithologic portion of the power spectra controlled by 9 datapoints for the longest line while the corresponding semivariogram has 131 datapoints of control.

The sparsity of control over the lithologic portion of the power spectra certainly is a drawback for utilizing this method for determining stochastic parameters. Even doubling the length of the longest line to 90 m would result in only one added datapoint. This contrasts with the well developed nested semivariograms seen on the semivariograms of 4 of the 5 seismic lines of this study. The only line difficult to interpret is line c in which the correlation length and the length of the line are nearly the same.

Observing the strong agreement of the stochastic parameters under study boosts confidence in use of the semivariogram for interpretation. Use of power spectra for interpreting these parameters is problematic. However,these plots can be used to confirm the interpretation obtained from the semivariograms. While the separation of the horizontal and vertical components of the power spectra help define the portions of the plots that can be utilized for lithologic interpretation, interpretation of the correlation length from these plots is difficult and merely enhances the results from the semivariogram method. While slopes corresponding to the fractal dimensions extracted from semivariograms and placed on the power spectra show reasonable agreement, close inspection of the power spectra show that a wide variation of slopes could have been used for this interpretation. Due to the scattering of the sparse data in the portion of the power spectra responding to lithology, a slope resulting in a fractal dimension of as low as 1.2 could have been used for this interpretation as well as those presented. This considerable variation in slope results in fractal dimension interpretations from the power spectra of questionable value other than for confirming the semovariogram method. Interpretation of power spectra have thus been utilized only to enhance the interpretation from the semivariogram methods.

While a local relationship between seismic reflectivity and hydraulic conductivity can be established, values of M for equation 9 should be established locally in areas of study. Local studies of hydraulic conductivity by collecting permeameter data from the face of study would contribute much to establishing this relationship. Establishing these values will contribute to predicting the distribution of hydraulic conductivities for the local area of study as well as the stochastic parameters of the formation under study.

The correlation length of a velocity field and its resulting hydraulic conductivity field calculated from equation 8 show good agreement (see Figs. 37 and 38). Further studies comparing correlation lengths of velocity fields and resulting hydraulic conductivity fields should be conducted to validate this relationship. While the difference between fractal dimension of the velocity field and the corresponding hydraulic conductivity field are not great, the differences of slopes seen on Figures 41 and 42 are enough to merit further study of the slope relationships. These studies could be accomplished utilizing known velocity fields for input.

Most semivariograms calculated from geologic field phenomena exhibit spherical model characteristics (this dataset being no exception after the lag distances greater than those corresponding to the wavelength of the seismic signal has been reached). Using this model on seismic attributes leads to knowledge of the percentage of variance from a particular point of interest within the field of interest. This could readily lead to contour maps of probability of similarity of the attribute of interest (or even a principal component of interest) around a particular point. If this point of interest were a well, the probability of an offset well exhibiting similar attribute(s) within given distances could be given for any point within the data volume. This application would be of interest in not only hydrogeologic studies, but also reservoir development studies in the petroleum industry.

Power spectra and semivariograms extracted from seismic data can be useful tools for interpreting the geostatistical parameters of correlation length and fractal dimension. These stochastic parameters can be related to the same stochastic parameters extracted from lithologic parameters of interest such as porosity, and hydraulic conductivity. These predicted stochastic parameters can then be used to determine lithologic variation within a field scale area of interest. While power spectra and semivariograms differ in method of extracting frequency (or wavelength) and correlation information, these differences only serve to enhance the interpretation when both methods are utilized. The semivariogram gives a considerably better value for both correlation length and fractal dimension and the power spectra shows the frequency bandwidth responding to seismic signal and lithologic variations.

Previous applications of extracting seismic attributes have included interpretation of rock type, porosity and even fluid type. While interpretation of absolute permeability values extracted from seismic data have not been utilized on as broad a scale as other properties, the use of the relative distribution properties of correlation length and fractal dimension on a local field scale shows promise (Mela, 1997). Advancement of techniques to extract these stochastic parameters related to their distribution will be of great value for future hydrogeologic studies.


For full sized graphics of figures (in pdf format) that can be printed out on 8 1/2 by 11 click here: (Fig.1, Fig.2, Fig.3, Fig.4, Fig.5, Fig.6, Fig.7, Fig.8, Fig.9, Fig.10, Fig.11, Fig.12, Fig.13, Fig.14, Fig.15, Fig.16, Fig.17, Fig.18, Fig.19, Fig.20, Fig.21, Fig.22, Fig.23, Fig.24, Fig.25, Fig.26, Fig.27, Fig.28, Fig.29, Fig.30, Fig.31, Fig.32, Fig.33, Fig.34, Fig.35, Fig.36, Fig.37, Fig.38, Fig.39, Fig.40, Fig.41, Fig.42, Fig.43, Fig.44, Fig.45)

Bear, Jacob, 1979, Hydraulics of Groundwater. New York, N.Y.: McGraw-Hill Publishing Co.

Carr, James R. and Myers, Donald E., 1984, Application of the Theory of Regionalized Variables to the Spatial Analysis of Landsat Data, IEEE Proceedings of the PECORA Information Technologies for Remote Sensing Today and Tomorrow

Carr, James R., 1995, Numerical Analysis for the Geological Sciences. Englewood Cliffs, NJ: Prentice Hall

Carr, James R., 1997, personal communication

Claerbout, Jon F., 1992, Earth Soundings Analysis Processing versus Inversion. Boston, Mass.: Blackwell Scientific Publications

Dobrin, Milton B., 1960, Geophysical Prospecting - 2nd Ed. New York, N.Y.: McGraw-Hill Book Co.

Dobrin, Milton B. and Savit, Carl H., 1988, Introduction to Geophysical Prospecting - 4th Ed. New York, N.Y.: McGraw-Hill Book Co.

Gelhar, Lynn W., Guttjahr, Allan L., and Naff, Richard L., 1979, Stochastic Analysis of Macrodispersion in a Stratified Aquifer: Water Resources Research, 15, 6, pp. 1387-1397

Gelhar, Lynn W., 1993, Stochastic Subsurface Geology. Englewood Cliffs, N.J.: Prentice Hall

Gibson, B.S., and Levander, A.R., 1988, Modeling and Processing of Scattered Waves in Seismic Reflection Surveys, Geophysics, 53, 4 pp 466-478

Knapp, R.W. and Steeples, D.W., 1986a, High-resolution Common-depth-point Reflection Profiling: Instrumentation, Geophysics, 51, 2 pp 276-282

Knapp, R.W. and Steeples, D.W., 1986b, High-resolution Common-depth-point Reflection Profiling: Field Acquisition Parameter Design, Geophysics, 51, 2 pp 283-294

Levander, Alan R., Hobbs, R.W., Smith, S.K., England, R.W., Snyder, D.B., and Hollinger, K., 1994, The Crust as a Heterogeneous "Optical" Medium, or "Crocodiles in the Mist": Tectonophysics, 232, pp. 281-297

Li, Dachang and Lake, Larry W., 1994, A Moving Window Semivariance Estimator, Water Resources Research, 30, 5, pp. 1479-1489

Louie, J. N., Clayton, R.W., and LeBras, R.J., Three-Dimensional Imaging of Steeply Dipping Structure Near the San Andreas Fault, Parkfield, California, Geophysics, 53, 2 pp. 176-185

Lindseth, Roy O., 1979, Synthetic Sonic Logs - A process for Stratigraphic Interpretation: Geophysics, 44, 1, pp 3-26

Mela, Ken, 1997, Viability of Using Seismic Data to Predict Hydrogeological Parameters, EEGS Proceedings of the Symposium on the Applications of Geophysics to Engineering and Environmental Problems, pp 773-780

Moore, James G., and Archibold, N.L. 1969, Geology and Mineral Deposits of Lyon, Douglas, and Ormsby Counties, Nevada, Nevada Bureau of Mines Bulletin 75

Press, W.H., Flannery, B.P., Teukolsky, S.A., and Vetterling, W.T., 1986 Numerical Recipes: The Art of Scientific Programming. Cambridge, UK: Cambridge University Press

Pullammanappallil, S.K., Levander, A., and Larkin, Steven P., 1996, Estimation of Crustal Stochastic Parameters from Seismic Exploration Data: submitted to Journal of Geophysical Research - July 1996

Schlumberger, Inc., 1972a, Log Interpretation, Charts: Houston

Schlumberger, Inc., 1972b, Log Interpretation, Vol. 1 - Principles: Houston

Steeples, Don W., Knapp, Ralph W. and McElwee, Carl D., 1986, Seismic Reflection Investigations of Sinkholes Beneath Interstate 70 in Kansas, Geophysics 51, 2 pp 295-301

Wheatcraft, S.W. and Tyler, S.W., 1988 An Explanation of Scale-Dependent Dispersivity in Heterogeneous Aquifers Using Concepts of Fractal Geometry, Water Resources Research, 24, 4, pp 566-578

Wyllie, M.R.J. and Rose, W.D., 1950, Some Theoretical Considerations Related to the Quantitative Evaluation of the Physical Characteristics of Reservoir Rock from Electrical Log Data: Journal of Petroleum Technology: 2, pp 105-118

Zhan, Hongbin and Wheatcraft, Stephen W., 1996 Macrodispersivity Tensor for Nonreactive Solute Transport in Isotropic and Anisotropic Fractal Porous Media: Analytic Solutions, Water Resources Research, 32, 12 pp 3461-3474