Resource Potential of Deep-Water Hydrates Across the Gulf of Mexico: Part 2, Evaluating Hydrate Systems with 4C OBC Seismic Data
Published:December 01, 2009
- PDF LinkChapter PDF
Bob Hardage, Diana Sava, Paul Murray, Mike DeAngelo, 2009. "Resource Potential of Deep-Water Hydrates Across the Gulf of Mexico: Part 2, Evaluating Hydrate Systems with 4C OBC Seismic Data", Unconventional Energy Resources: Making the Unconventional Conventional, Tim Carr, Tony D’Agostino, William Ambrose, Jack Pashin, Norman C. Rosen
Download citation file:
We have evaluated hydrate concentrations across deep-water areas of Green Canyon, Gulf of Mexico, using well log data and four-component (4C) seismic data acquired by companies interested in deep oil and gas targets, not in near-seafloor hydrates. Even though the data are not acquired for purposes of studying nearseafloor geology, we have found these off-the-shelf industry data to be invaluable for evaluating hydrate systems positioned close to the sea floor. We summarize our data analyses and initial research findings in this publication as a two-paper sequence.
In this second paper of our two-part series, we describe how two images of deep-water hydrate systems can be made from four-component ocean-bottomcable (4C OBC) seismic data: a compressional (P-P) image and a converted-shear (P-SV) image. We further illustrate how we implement a raytracing procedure to determine accurate values of P-wave velocity (VP) and SV-mode velocity (VS) across thin subsea-floor layers. These interval velocities are used with the rock physics theory described in our first paper of this two-paper sequence to (1) estimate hydrate concentration at calibration wells where there are resistivity logs to use for an independent calculation of the hydrate fraction, and (2) expand the hydrate estimation along 4C seismic profiles that extend long distances away from calibration wells.
We present maps of hydrate concentration estimated by our data analyses and physical assumptions. We found hydrate concentration to not exceed 40 percent of the available pore space across our Green Canyon study area and to usually be in the range of 10 to 20 percent of the available pore volume.
We developed a new seismic data-processing concept that creates the down-going compressional (P) wave field and the upgoing P and converted-shear (SV) wavefields that are needed to construct P-P and P-SV images of near sea-floor geology. In our data-processing approach, we worked with four-component ocean-bottom-cable (4C OBC) data collected at each individual sea-floor receiver station, segregated these data into common-receiver gathers for each of the four sensor components, and then used these four common-receiver gathers to create the downgoing P wavefield and the upgoing P and SV wave fields at that sea-floor receiver location. With these three wavefields (P-down, P-up, and SV-up), we then proceeded to calculate P-P and P-SV reflectivities of near-seafloor strata (Backus et al., 2006; Hardage et al., 2009).
Our data-processing concepts can be applied to 4C OBC data acquired in deep water only. Our data-processing assumptions are invalid when the vertical distance between the source and receiver (i.e. water depth) is of the order of 200 m or less. There is a strong similarity between our processing of sea-floor based common-receiver gathers and the processing of walkaway vertical seismic survey (VSP) data. No one in the seismic data-processing industry has utilized this deep-water 4C OBC data-processing technique before to our knowledge.
4-C Seismic Data
In our seismic imaging approach, 4C OBC seismic data was treated as common-receiver gathers. The 4C data acquired at one typical sea-floor receiver station inside our study area are shown as Figure 1. All of our data processing was performed on individual receiver-station data like the data illustrated in this figure. Receiver stations were spaced at 25-m intervals along all of the OBC profiles that traversed our study area, and the data-processing procedure described here was repeated at each of these receiver stations.
The fundamental theory of our P-P data-processing strategy is based on analyzing data that have been acquired with a sensor that has a hydrophone and a vertical geophone. The key sensor-response equations involved in OBC data acquisition are illustrated and explained on Figure 2.
Defining D as the downgoing compressional wavefield that reaches a seafloor station and U as the upgoing compressional wavefield at that same station, the responses of the hydrophone P, vertical geophone Z, and inline horizontal geophone X are:
Φ is the incident angle at which the downgoing compressional wave arrives at the seafloor (Fig. 2). These equations imply that after appropriate calibration, a seafloor hydrophone response (P) and a seafloor vertical-geophone response (Z) can be combined to create the unknown down-going (D) and up-going (U) P-P wavefields at each receiver station using the following relationships:
By having access to the down-going (D) and up-going (U) compressional wave fields, subsea-floor P-P reflectivity R can be recovered by taking the ratio,
in the frequency domain, f. The inverse Fourier transform of RPP(f) then creates a time-based reflectivity series that starts at the seafloor and extends to a depth below the base of the hydrate stability zone. It is this time-based reflectivity RPP(t) that we use to create our high-resolution images of deep-water, near-seafloor geology.
We first calculate the Fourier transform of the up-and down-traveling P-P waves obtained from the simple combinations of raw P and Z data that are defined by Equations 1 through 5. At each offset, we divide the up-traveling P-P wave by the down-traveling P wave, using a modest damping applied for stability. An inverse Fourier transform then yields the P-P reflectivity. The reduction to sea-floor datum is automatic in this ratio process.
To determine the reflectivity of the P-SV wave, we follow a similar procedure. We divide the extracted P-SV wave by the down-going P wave field (lower left, Fig. 3) in the frequency domain to produce P-SV reflectivity defined as:
An inverse Fourier transform produces the P-SV reflectivity. More detail about this technique for calculating P-P and P-SV reflectivities has been published by Backus et al. (2006) and by Hardage et al. (2009).
Creating Local Common-Receiver Images
To create P-P and P-SV images, we apply dynamic corrections to our reflectivity estimates to correct for the moveout on the near-offset traces. On Figure 3, we show the deconvolved P-P data (that is, the P-P reflectivity calculated by the strategy that leads to Eq. 6) for a ±2,500-m offset range after applying a time differentiation to enhance the frequency of the data. The data are excellent quality over the full offset range. Raytracing using a layered-velocity model of subsea-floor geology is then used to calculate curves of source-receiver offset versus time that correspond to reflection depth points that are a fixed offset distance from the receiver location. Examples of raytrace curves calculated for this common-receiver gather are shown in Figure 4 for depth-point offsets starting at ±10 m from the receiver station and increasing at 25-m intervals out to ±160 m from the receiver coordinates. Image-trace data can now be recovered by interpolation along these curves to produce P-P image traces at specified depth-point offsets from the receiver location.
In Figure 5, we show the deconvolved P-SV common-receiver gather at this same receiver station, before and after dynamic corrections. Note that even before dynamic moveout corrections, the P-SV events are nearly flat, so a limited-range stack before applying a dynamic moveout correction can provide a fairly good P-SV image. Depth-point-offset curves overlain on the P-SV reflectivity show that for our OBC data, which have a high VP/VS velocity ratio, any P-SV image trace of near-seafloor geology will extend only 1 or 2 m away from a receiver station.
Creating Continuous Images Along an OBC Profile
When these data-processing steps are followed at all receiver stations along an OBC profile, miniscale P-P and P-SV images are created at each receiver station. Each mini-image we created represents the subsea-floor image across a 25-m distance centered on a receiver station, which is the receiver-station interval for the OBC data that we used in this study. We then combine these small-scale images to make continuous P-P and P-SV images that extend for several kilometers along each OBC profile. This concept is illustrated on Figure 6 for P-P imaging and as Figure 7 for P-SV imaging.
The 4-km P-P image shown on the left of Figure 6 is a series of small-scale, local P-P images constructed at each sea-floor receiver station. Each local image is created by first calculating P-P reflectivity at each receiver station using Equation 6. The P-P reflectivity at one arbitrary sea-floor station A is shown on the right of the figure. We then determine a sequence of constant-depth-point-offset traces across these P-P reflectivity data, such as a data trace along any of the red curves that overlay the P-P reflectivity. We arbitrarily decided to interpolate a P-P image trace at depth-point-offset intervals of 5 m. Five of these traces create a 25-m-wide P-P image centered on the receiver station. For the data on Figure 6, we created a five-trace, 25-m-wide, local image at 160 consecutive receiver stations to make the 4-km P-P image that is displayed on the left.
This same concept is repeated in Figure 7 for the P-SV reflectivity. The principal difference between P-SV imaging and P-P imaging is that the low VS velocities in deep-water, near sea-floor strata do not allow constant-offset depth-point image traces to be calculated at distances farther than 1 or 2 m from each seafloor receiver station. As a result, we could not create a 25-m-wide P-SV image using 5-m trace spacing as we did with P-P data. Instead, we create a single, zero-offset P-SV image trace at each receiver coordinate by summing all the traces between the +1 m and -1 m depth-point-offset curves (Fig. 7). The result is a P-SV image along each OBC profile that had a trace spacing of 25 m, the same distance as the receiver-to-receiver interval.
Comparison with State-of-the-Art Imaging
On Figure 8 we compare a P-P section produced by state-of-the-art imaging done by a leading seismic contractor (left) with our P-P imaging technique that utilizes common-receiver gathers (right). The data are displayed with a seafloor datum. The ocean-floor multiple arrives at a time greater than the maximum time shown on this P-P section, which allows us to use our data-processing simplifications to produce this excellent-quality image. The improvement in detail in our image compared with that of conventional processing is striking. In addition to improved vertical resolution, there is a marked increase in structural detail and horizontal resolution, even though the contractor data have been migrated and our data have not. Of special interest is the comparison on the far left of each image, especially above 100 ms, where our image shows strata that dip sharply to the left, the same stratigraphic dip observed in the kilohertz-range P-P data that are exhibited later.
On Figure 9 we compare a state-of-the-art P-SV section (produced by the same leading contractor that made the P-P image in Figure 8) with our processing approach that is based on common-receiver gathers. We see that our technique provides significant data quality improvement close to the sea floor. Note the strong P-SV reflection at a shallow depth of 1.5 m that parallels the sea floor. The strong reflection at a subsea-floor depth of 10 m is at the center of a sedimentary package that is unconformable with the sea floor.
What we consider amazing is the resolution of the P-SV images produced by our imaging technique. To illustrate this resolution, we compare on Figure 10 our air-gun-frequency (10‒200 Hz) P-SV image with a high-frequency (2‒10 kHz) P-P image acquired by deep-water Autonomous Underwater Vehicle (AUV) technology. Depth-equivalent P-P and P-SV horizons are identified on the images. These images show that the frequency of P-P data has to be increased into the kilohertz range for P-P data to have a spatial resolution that is equal to that of low-frequency P-SV data. This observation confirms the simple principle that in order for P-P data and P-SV data to have equivalent resolution, they must have equivalent spatial wavelengths. This requirement leads to the simple conclusion that if
where λ is wavelength and subscripts p and s refer, respectively, to P-wave and S-wave, then the frequencies (f) must be related as,
This wave physics requires P-P frequency to be boosted by an amount equal to VP/VS relative to P-SV frequency in order for P-P data to have the same wavelength spectrum as P-SV data. In the near sea-floor strata across our study area, VP/VS can exceed 55 (Fig. 10, top layer), which requires that P-P frequency be boosted into the kilohertz range to match the resolution that we achieve with 10 to 200 Hz P-SV data using our data-processing strategy.
Interpretation of OBC Profiles
The objective of our interpretation of P-P and P-SV images along each OBC line was to define which subseafloor P-SV reflection events were depth-equivalent to selected P-P reflections and to calculate VP/VS velocity ratios within the stratigraphic intervals bounded by these depth-equivalent P-P and P-SV horizons. We limited our interpretation to reflection events that extended only a short distance below the base of the hydrate stability zone (BHSZ), such as the horizons shown on the profiles exhibited as Figure 11. We estimated the subseafloor depth of the BHSZ from the prediction guidelines developed by Milkov and Sassen (2001) for the Green Canyon area. Once the depth position of the BHSZ was estimated from these Milkov and Sassen calibration curves, we then used logical velocity assumptions to interpret a pair of depth-equivalent P-P and P-SV horizons positioned below that estimated depth. These horizons were the deepest interfaces that we used in our seismic velocity analysis.
We concluded that determining depth-equivalent subseafloor units and making reliable velocity estimates within near-seafloor layers required that a rigorous numerical raytracing analysis be done to determine if each pair of tentatively interpreted P-P and P-SV reflections were truly depth-equivalent, or whether different events needed to be selected to establish depth equivalency. The following section illustrates how this raytracing procedure was performed at the receiver station located where the vertical line extends across the P-P and P-SV images displayed on Figure 11.
The concept of the raytracing method we implemented is illustrated by the diagram on Figure 12. At intervals of 10 sea-floor receiver stations (250 m) along each OBC profile, a common-receiver gather was constructed using data generated at 120 source stations centered on the receiver position (-3000 m to +3000 m source offsets). When the P-P and P-SV images showed geology along a profile to be laterally uniform, the distance between consecutive raytrace analysis points was sometimes extended to span 20 receiver stations (500 m).
Once depth-equivalent horizons are interpreted on stacked P-P and P-SV data, as has been done on the P-P and P-SV images along the profile on Figure 11, VP/VS velocity ratios can be calculated between any two horizon pairs using the equation,
where Δ TPS is the interval time between the two P-SV horizons that bound the interval on the stacked P-SV image, and Δ TPP is the interval time between the two corresponding P-P depth-equivalent horizons on the stacked P-P image. The accuracies of the VP/VS values determined in this manner depend on the skill of the interpreter in recognizing depth-equivalent P-P and P-SV reflectivity behavior.
Our raytracing methodology uses these interpreter-defined VP/VS velocity ratios as a fixed Earth property. As the raytracing progresses from the first subseafloor layer to the last layer that extends below the base of the hydrate stability zone, the philosophy of our raytracing approach to building an Earth model of subseafloor velocity layers is: (1) assume the VP/VS ratios determined by interpreting depth-equivalent P-P and P-SV reflections across each subsea floor interval are correct for all intervals, (2) use each raytrace-based interval value of VP (or VS) together with the fixed VP/VS velocity ratio across the interval to determine the corresponding interval value of VS (or VP), and (3) define the picked reflection time of the imaged horizon at the receiver-station coordinates where the velocity analysis is done as the zero-offset time that must be associated with the proper reflection event in the common-receiver gather. Using these concepts, we raytrace down from each sea-level source to the target horizon and then up from that horizon to the sea-floor receiver and adjust the velocity model until raytrace-calculated arrival times converge to the observed reflection arrival times of the proper event in the P-P and P-SV common-receiver gathers. The concept is depicted graphically on Figure 13.
Because raytracing is not computationally intensive, we can iterate the arrival-time calculation using different layer thicknesses and layer velocities and perform each velocity analysis interactively. The option shown on Figure 13 of reinterpreting the stack data to refine the definition of depth-equivalent P-P and P-SV reflections is important for accurate velocity analysis because the parameters T0(PP), T0(PS), and VP/VS associated with these depth-equivalent events are the critical constraints used in the raytracing analysis. Stack data are sometimes reinterpreted three or four times to achieve raytrace convergence of VP and VS across a layer at some receiver stations.
In practice, we can overlay raytrace-based travel-time curves on the common-receiver gathers, or we can apply raytrace-based time shifts to each trace of the gather to flatten a targeted reflection that is being analyzed.
We found the latter approach (time shifting and reflection flattening) to be simpler and more accurate. It is important to note that flattening an event by static shifts is not equivalent to applying dynamic moveout corrections to the data. The flattening process simply shifts each data trace by a calculated time, thereby eliminating wavelet-stretch artifacts and allowing us to use much longer source-receiver offsets to constrain velocities than can be done in a time-variant, moveout-based approach. An additional benefit of a raytracing approach is that it is not susceptible to limitations of nonhyperbolic moveout associated with low velocities, large offset/depth ratios, and differences in source and receiver elevations. All three of these conditions are present in the deep-water OBC data used in this study.
A principal advantage of our velocity analysis strategy is that velocity layers are defined as a function of depth below the seafloor. As a result, VP and VS velocities derived by an analysis can be correlated with depth-based resistivity logs, and both resistivity and velocity data can be used to identify the BHSZ boundary in P-P and P-SV image spaces. This Earth-layer construction process has been applied at approximately 800 seafloor receiver stations across the grid of 2D OBC lines to build a continuous velocity model along each 2D OBC profile. Velocity layer 1 starts at the sea-floor and extends to the shallowest interpretable P-P reflection. Velocity layers 2, 3, and 4 extend to successively greater depths until a velocity layer N is created that extends deeper than the BHSZ boundary.
Integration of Resistivity, Velocity, and Seismic Data
The resistivity and velocity profiles at one calibration well will be used to demonstrate how these Earth properties can be correlated with P-P and P-SV images and seismic attributes along OBC profile 264. First, the layer-velocity model built at this well is adjusted to match the P-P and P-SV image-time axes at the well location, as shown in Figure 14. This correlation process allows depth-based data to be compared against time-based seismic velocities. Inspection of the figure shows that each Earth-velocity layer correlates with a distinct seismic facies unit in both P-P image space and in P-SV image space.
Three estimates of the base of the hydrate stability zone [labeled BHSZ(90%), BHSZ(R), BHSZ(V)] are marked on the seismic profile. These horizons have the following meanings:
BHSZ(90%): The depth of the base of the hydrate stability zone for a natural gas having 90.4 percent methane as calculated by Milkov and Sassen (2001).
BHSZ(R): The depth of a decrease in formation resistivity that is “close to” the depth of horizon BHSZ(90%) and that appears to be a logical choice for the base of the hydrate stability zone when examining resistivity log data acquired in the calibration well.
BHSZ(V): The depth of a decrease in VP velocity that is “close to” the depth of horizon BHSZ(90%) and that appears to be a logical choice for the onset of free-gas trapped below the base of stable hydrate, as defined by seismic VP interval velocities.
It is important to note that the VP velocity profile at this well exhibits an increasing trend in magnitude through Layer 4 and then undergoes a velocity reversal in Layer 5. A tentative dilemma presented by this data-correlation exercise is that formation resistivity increases in Layer 5, indicating increased hydrate content in that layer; whereas, the P-wave velocity decreases, which indicates decreased (or absent) hydrate content. We thus have opposing interpretations: resistivity data imply hydrate is present in Layer 5, but velocity data indicate hydrate is absent.
We interpret the increased formation resistivity in velocity Layer 5 to be caused by free gas, not by hydrate. This adjusted interpretation of the resistivity log brings the resistivity data and velocity data at this calibration well into agreement because the decrease in VP velocity in Layer 5 is also consistent with the presence of free gas. From this logic, we adjust the base of hydrate stability upward to depth BHSZ(V), the base of velocity Layer 4 where the reversal in VP velocity begins.
2D Profiles of Velocity Layering
After performing joint inversions of log-based resistivity and seismic-based velocity at several calibration wells (paper 1 of this two-part series), we determined an optimal function that could be used to relate hydrate concentration to seismic-based VP velocity at OBC line coordinates between calibration wells. The input data for this velocity-based hydrate estimation were 2D profiles of VP layer velocities determined by raytrace analysis of common-receiver gathers (see Fig. 13). These raytrace analyses were done at intervals of 10 receiver stations (250 m) along each OBC profile. Because there were approximately 8,000 receiver stations across the OBC grid we studied, our approach to velocity analysis required that we construct more than 800 models of depth-based layers of VP and VS velocities along the approximately 200 km of OBC profiles that we analyzed. This velocity study was laborious and time consuming, but was essential for reliable hydrate estimation. An example of the types of 2D velocity layer models we created is exhibited as Figures 15a and 15b.
We made no attempt to smooth the velocity values displayed on this figure (or on any of the other velocity profiles) in order to make data displays have a more pleasing cosmetic appearance. Our velocity analyses across the OBC data grid consistently indicated a VP velocity inversion occurred near the depth of the expected BHSZ boundary, as illustrated on Figure 14.
2D Profiles of Hydrate Concentration
Relationships between VP velocity and hydrate concentration developed at calibration wells were applied to the VP velocity layer models constructed along each OBC profile. The inversion results for the velocity layering along profiles 549 and 553 are displayed as Figure 16. Along the southern half of each profile, the BHSZ boundary was defined as the onset of a reversal in VP magnitude. Along the northern half of each line, the BHSZ boundary was defined by the water-depth-based thermal constraint for 90 percent methane hydrate published by Milkov and Sassen (2001).
Because the normal compaction curve has such a dynamic depth variation across velocity Layer 1 immediately below the seafloor and no log data is available across this shallowest layer to confirm the effect of compaction on velocity, we have assigned a constant, near-zero hydrate concentration to Layer 1 and focus our hydrate estimation on velocity Layer 2 and deeper layers that extend down to the BHSZ horizon. Our velocity analyses does not indicate a velocity magnitude in Layer 1 anywhere across the OBC profile grid that would infer hydrate is present in this shallowest layer.
Our calculated hydrate concentrations exhibit considerable lateral spatial variation within each velocity layer and even greater vertical variability from layer to layer. The maximum hydrate concentration found along the two particular OBC profiles exhibited as Figure 16 were local, limited areas where hydrate occupied a little more than 30 percent of the pore space of the host sediment.
Mapping the Amount of In Situ Hydrate
To determine the amount of in situ hydrate existing within the interval extending from the seafloor to the BHSZ boundary, we multiply our seismic-based hydrate concentrations (expressed as the fraction of occupied pore space) by each layer thickness and layer porosity and sum these products to create an estimate of total in-place hydrate. The resulting maps of in-place hydrate across the study areas are shown as Figures 17 and 18. A regional map showing the locations of these two study areas is provided in Figure 1 in Part 1 of this two-part series (Sava and Hardage, this issue). Our seismic-based quantification of in situ hydrate indicates the largest accumulation of hydrate within our limited study area exists north of Genesis Field (Fig. 18).
At some locations across this trend, the amount of in-place hydrate is estimated to be as much as 2000 to 4000 m3 beneath 1 × 250-m rectangular strips centered on receiver stations where VP interval velocities were determined for estimating hydrate concentration. Other significant accumulations of hydrate are shown by the green to red colors that are shown at several locations across the OBC grid.
Comparing Load-Bearing and Free-Floating Hydrate Assumptions
The hydrate distributions displayed as Figures 17 and 18 are estimated using the assumption that the hydrate granules embedded in the sediment bear a proportionate part of the sediment weight. This assumption leads to the “load-bearing” rock physics theory of hydrate morphology analyzed in paper 1 (Sava and Hardage, this issue) of this two-paper series. An alternate assumption that has merit is that unit volumes of hydrate float in the pore spaces of the host sediment and are not part of the load-bearing matrix. This assumption leads to the “free-floating” rock physics theory of hydrate morphology, also discussed in paper 1.
For a given value of VP within a near-seafloor layer, a free-floating assumption for the hydrate-sediment morphology results in greater hydrate saturation than does a load-bearing assumption. A comparison of the hydrate concentrations predicted by these two hydrate-morphology models along OBC profile 557 is displayed as Figure 19. For the range of interval VP velocities that we found within the hydrate stability zone in the Green Canyon area, our free-floating-hydrate theory causes approximately five more percentage points to be added to the hydrate fraction than what is predicted by our load-bearing-hydrate theory. The almost-constant difference of approximately five percentage points of hydrate concentration that results when using these two hydrate-morphology assumptions is illustrated by the profiles displayed as Figure 19.
We emphasize that our seismic data-processing methodology is specialized for optimal imaging of deep-water, near sea-floor geology. Because our imaging procedure is to treat a seafloor-receiver trace gather as a walkaway vertical seismic profile, there needs to be a significant vertical distance between the surface-based source and the seafloor receiver. We have not applied our data-processing strategy to 4C OBC data acquired in water depths less than 450 m and do not know the minimum water depth where our data-processing technique begins to produce unacceptable results.
We do not champion the use of our data-processing technique for imaging geology at significant depths (≥1000 m) below the sea floor. Although our data-processing concepts should produce acceptable quality images of geology at subseafloor depths of 1000 m and more, our procedure appears to have no advantage for imaging deep geology over that provided by the procedures used across the seismic data-processing industry. The decades-old common-midpoint (CMP) procedures used in the industry are superb for creating P-P images of deep geology. The common-conversion-point (CCP) binning/imaging techniques used by commercial seismic data-processing shops for generating P-SV images is also an excellent method for imaging deep geology. Our imaging strategy is specialized to produce optimal images of near sea-floor geology where deep-water hydrate systems occur.
One limitation of our current data-processing method is that we make no attempt to remove water-column multiples from the seismic data. Once these multiples begin to appear in the reflection data, there is often a reduction in the signal:noise ratio across some depth intervals and less-accurate imaging of the subseafloor geology. Because the water depth ranged from 500 m to 1000 m along the OBC profiles that we analyzed, water-column multiples do not appear in the data until the P-P image extends below the base of the hydrate stability zone. We believe we attenuated water-column multiples adequately for the deeper parts of our P-SV images where some multiple contamination did exist. We can expand our data-processing procedure to attenuate water-column multiples if 4C OBC data analyses are done at other deep-water sites, but we do not see the need to do so for the Green Canyon area that we studied. Published papers supporting or complementing the research findings reported herein would be studies by Carcione and Tinivella (2000), Ecker and others (2000), Lu and McMechan (2002), Miller and others (1991), Wood and others (1994), and Zillmer (2006).