Patent Publication Number: US-9405026-B2

Title: Estimation of production sweep efficiency utilizing geophysical data

Description:
CROSS-REFERENCE TO RELATED APPLICATION 
     This application claims the benefit of U.S. Provisional Patent Application 61/569,576 filed Dec. 12, 2011, entitled ESTIMATION OF PRODUCTION SWEEP EFFICIENCY UTILIZING GEOPHYSICAL DATA, the entirety of which is incorporated by reference herein. 
     FIELD OF THE INVENTION 
     This invention relates generally to the field of oil and gas production and, more particularly, to estimating production-induced changes in hydrocarbon fluid saturation and the spatial distribution of the residual hydrocarbons utilizing geophysical information. Specifically, the invention is a method for characterizing saturation heterogeneities in reservoirs (i.e., unproduced pockets of oil or gas) by quantifying and interpreting changes in P-wave velocity and attenuation utilizing rock physics relationships. 
     BACKGROUND OF THE INVENTION 
     The objective of the present invention is to use attenuation and velocity dispersion data to determine parameters that describe subsurface fluid distribution. These parameters may include spatial distribution fluids, heterogeneity length scale, and average saturation. Specifically, the invention employs a poroelastic rock physics model (and its frequency behavior or asymptotes) to provide a physics-based interpretation of geophysical parameters. Published work that might be relevant to the invention falls into two classifications. They are either theoretical studies or patents which seek to invert geophysical data into fluid distribution parameters, mainly average saturation. See U.S. Pat. No. 7,136,757 to Goloshubin and Korneev; U.S. Patent Application Publication No. 2008/0175099 by Hawthorne et al.; and Tserkovanyak and Johnson, 2004). There are many theoretical modeling studies that demonstrate aspects of the physical mechanism employed by the invention (see, for example, Johnson 2001; Pride at al. 2004; Toms et al. 2008), but none of these studies present any strategy to interpret the geophysical information for parameters describing the fluid distribution. The present inventors are not aware of any published work aimed at using attenuation to define all three significant aspects of subsurface fluid distribution (spatial distribution, length scale and average saturation). 
     In U.S. Pat. No. 7,136,757, Goloshubin and Korneev (2006) utilize the reflection characteristics of seismic waves from a target layer and reference layer of known average fluid saturation to determine the average fluid saturation within the target layer. The reference layer is normally decided from the borehole or log data; however other reference layers may be utilized. They perform a low pass filtering of the windowed reflections from the target layer leaving frequencies below the lowest recorded frequency spectra. They define a complex ratio R(x,ω)=W(x,ω)/W(x0,ω) of the low frequency amplitude spectra of the target W(x,ω) and reference layer W(x0,ω) data. They determine the average relative reflected amplitude A(x) by summing the ratios R(x,ω) of the first N reliable low frequency ω. They also analyze the derivative of the average amplitude D A (x) with respect to frequency and the derivative of the average time delay D P (x) with respect to frequency. They map the average fluid saturations by utilizing the calibration conditions: A(x 0 )=1, D A (x 0 )=0, and D P (x 0 )=0 at the location x=x 0 . The invention patented by Goloshubin and Korneev estimates only average fluid saturation and fluid type based on attribute (seismic spectral) variation from the calibration conditions. While attenuation and phase velocity dispersion may be implicit in the physics of this invention, in that the reflected amplitude from the target and reference layers will be affected similarly by velocity dispersion, this additional physical knowledge is not exploited to determine the full gamut of subsurface fluid distribution parameters. 
     In U.S. Patent Application Publication US 2008/0175099, Hawthorne et al. determine fluid type within subsurface formations surrounding a borehole. The invention utilizes sonic logging data, in addition to open-hole log measurements of porosity, lithology etc. They propose three different ways of characterizing fluid type. The first is to determine formation mobility, which is the ratio of rock permeability to fluid viscosity. They propose using the formation mobility to indicate fluid viscosity which differs for heavy oil, oil, gas and water. The second approach considers pore fluid bulk modulus, as the compressibility of heavy oil, water, gas etc is very different. The third approach considers attenuation. The invention utilizes compressional, shear and Stoneley waveforms, where slowness (the inverse of velocity) and attenuation of each waveform may be utilized to invert for pore fluid bulk modulus, borehole mud slowness or formation mobility. 
     Tserkovank and Johnson (2001) seek to invert geophysical data for heterogeneity length scale information. They utilize the laboratory experimental data acquired by Cadoret et al. (1998). They assumed a respective water saturation, and for each saturation used the corresponding P-wave attenuation and velocity measurements to determine a length scale for the heterogeneity. Their method employs a rock physics model that assumes a fixed fluid distribution pattern described by periodically distributed spherical saturation heterogeneities. The approach described in the Tserkovanayak and Johnson paper employs a straight inversion of attenuation and velocity information to determine length scale, utilizing a rock physics model knowing all other rock and fluid properties, including the average fluid saturation. 
     SUMMARY OF THE INVENTION 
     In one embodiment, the invention is computer implemented method for estimating fluid heterogeneity in a subsurface region from compressional wave attenuation or velocity dispersion, comprising: (a) measuring compressional wave attenuation or velocity, or extracting it from geophysical data, for at least one frequency; (b) selecting a frequency-dependent, mathematical rock physics model for attenuation or velocity, said model pertaining to or containing one or more model parameters that pertain to heterogeneous features of fluid distribution in the subsurface region; (c) using the rock physics model to predict frequency dependence of compressional wave attenuation or velocity; and (d) using the predicted frequency dependence, the attenuation or velocity measured or determined from the geophysical data for at least one frequency, and measured or estimated rock and fluid properties of the subsurface region, called subsurface properties, to determine, using a computer, at least one unknown subsurface property, including fluid saturation, that is related to fluid heterogeneity, each such subsurface property corresponding to a parameter in the mathematical rock physics model. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which: 
         FIG. 1  shows two examples of potential hydrocarbon spatial distribution within an otherwise water-saturated rock; 
         FIG. 2  is a flowchart showing basic steps in one embodiment of the invention; 
         FIGS. 3A-3C  show P-wave velocity dispersion for rocks of different permeability; 
         FIGS. 4A-4C  show frequency dependent P-wave attenuation for rocks of different permeability; 
         FIGS. 5A-5C  show P-wave dispersion for rock of permeability 1 Darcy, for varying length scales of saturation homogeneities; 
         FIGS. 6A-6C  show frequency dependent P-wave attenuation for rock of permeability 1 Darcy, for varying length scales of saturation homogeneities; 
         FIG. 7  shows peak frequency of attenuation for different permeability rocks having different heterogeneity length scales, where frequency is represented by gray-scale shading; 
         FIG. 8  is an attenuation cross plot of fluid contrast and rock compressibility, where the gray-scale is a measure of attenuation; 
         FIG. 9  is a cross plot like  FIG. 8 , except that the ratio of concentration of the most compressible fluid to the least compressible fluid is 0.4 instead of 0.1 as it is for  FIG. 8 ; 
         FIGS. 10A-10B  show attenuation (A) and velocity dispersion (B) estimated using a Mueller and Gurevich rock physics model by assuming two different spatial distributions of fluids, random and periodic; 
         FIG. 11  is a plot of attenuation versus frequency showing a high-frequency asymptote and a low-frequency asymptote derived from the Mueller and Gurevich rock physics model according to the present invention; 
         FIG. 12  is a plot like  FIG. 11  showing key features of utilizing the asymptotes to interpret attenuation behavior according to the present invention; 
         FIG. 13  shows the inverted asymptotes plotted as a cross plot between heterogeneity length scale and percentage of water saturation, using the data on attenuation versus frequency that was modeled for this example, where the intersection of the two inverted asymptote defines the length scale and relative saturation of the fluid distribution; and 
       for  FIGS. 14-16 , the heterogeneity scale length as a function of water saturation from  FIG. 13  is substituted into the rock physics model, and for the frequency corresponding to each of the three input data points, attenuation is plotted versus percentage of water saturation for both the low-frequency asymptote and the high-frequency asymptote of the model, enabling a determination of whether each input data frequency point is located at the peak frequency of attenuation or at frequencies above or below the peak frequency. 
     
    
    
     The invention will be described in connection with example embodiments. However, to the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claim. 
     DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS 
     This invention relates to estimating original or production induced changes in hydrocarbon fluid saturation and the spatial distribution of the residual hydrocarbons utilizing geophysical information. These production measurements may be interpretable by quantifying changes in P-wave velocity and attenuation utilizing rock physics relationships. This invention aims at characterizing saturation heterogeneities in reservoirs that are significantly smaller than the P-wave wavelength, which cannot normally be resolved utilizing standard geophysical methods. Assuming a rock physics model, the expected behavior of attenuation and/or velocity as a function of frequency can be predicted. The model-based predicted behavior along with attenuation and/or velocity measurements at one or more frequencies and assumptions of the reservoir conditions (temperature, pressure) and properties (rock and fluid type, elastic properties, etc.) are used to constrain the type of spatial distribution, length scale of the heterogeneities, or fluid saturation. 
     The results of the invention can be utilized for refining initial models of reservoir fluid saturation prior to staged production simulation and for refining estimates of hydrocarbon volumes still in place. In the production setting, the results of the invention may also be utilized to diagnose sweep efficiency. That is, whether fluid sweeping across the reservoir is homogeneous or heterogeneous and if heterogeneous at what length scale. The results of the invention may also be used in an exploration setting for hydrocarbon prediction; however there is more uncertainty in this scenario, due to the lack of prior rock property knowledge. 
     When producing from reservoir rock exhibiting spatial variation in permeability, complexity in wettability or spatial variation in pore topology, a significant volume of hydrocarbon may remain in place. To assist with maximizing recovery, it is beneficial to monitor production induced changes in fluid saturation utilizing geophysical data. The usual approach is to conduct a geophysical survey both prior to production (base survey) and subsequent to production (monitor survey). When survey and data processing repeatability is high, the differences between geophysical data sets may be attributed to alterations in reservoir pressure or fluid saturation. As of yet, such time-lapse, sometimes referred to as 4D, signatures observed in geophysical data are not readily interpretable in a quantitative way or even in a qualitative way for certain reservoir rock types. In part this may be a consequence of not taking into account additional physical effects, which result from the presence of sub-wavelength fluid saturation heterogeneities. Here saturation heterogeneities refer to immiscible (interfacially separated) fluids, which are spatially distributed through out the rock non-uniformly. Sub-wavelength refers to a characteristic length scale of the spatial distribution being less than the wavelength of the seismic compressional wave. The presence of sub-wavelength saturation heterogeneities has been shown experimentally and theoretically to cause P-wave attenuation and velocity dispersion (velocity dispersion refers to the variation of wave speeds with frequency). See, among other references, Cadoret et al. (1995), Cadoret et al. (1998), Johnson (2001), and Müller and Gurevich (2004). 
     Common Fluid Substitution Practices Utilized in 4D Studies 
     Crucial for interpreting 4D signatures for changes in fluid saturation is the capacity to model the elastic properties of partially saturated rock. (Partial saturation here means that the rock has at least two pore fluids, typically oil and brine.) Most techniques for interpreting or forward modeling the elastic properties of the fluid saturated rock assume a uniform distribution of fluids and utilize Gassmann&#39;s fluid substitution relationship (Gassmann 1951). By knowing the elastic properties of the pore fluid, dry rock and grain that makes up the rock, in addition to the porosity of the rock, Gassmann&#39;s relationship can be utilized to calculate the elastic properties of a fully fluid saturated rock. The influence of multiple pore fluids is then addressed by substituting an effective pore fluid with an elastic modulus that is the volumetric average of the elastic moduli of the component pore fluids. 
     There are three effective pore fluids commonly substituted into Gassmann&#39;s relationship for evaluation of 4D signatures. The first approach assumes that pore fluids are uniformly distributed through out the pore space, or that within a half wavelength of the P wave, that fluid pressure is equilibrated. When either of these assumptions is satisfied, the bulk modulus of the effective pore fluid is equal to the harmonic volume average of the elastic moduli of each pore fluid. This effective pore fluid is then utilized within Gassmann&#39;s relationship to model the bulk and shear moduli. 
     The second approach assumes that pore fluid does not equilibrate within a half wave cycle. This is modeled by assuming that the shear modulus of the rock is independent of the elastic properties of the pore fluids. When this assumption is permitted, the bulk modulus of each region of the rock saturated by a different pore fluid is modeled using Gassmann&#39;s fluid substitution method. The overall elastic properties of the rock are then modeled assuming Hill&#39;s relationship (Hill 1963), which relates the inverse saturated P wave modulus to the harmonic volume average of the saturated P wave modulus of each rock component saturated by a different pore fluid. Commonly the second approach is approximated by utilizing Gassmann&#39;s fluid substitution and postulating the presence of an effective pore fluid, which has a bulk modulus equivalent to the arithmetic volume average of the bulk modulus of each pore fluid. 
     The third approach also postulates the existence of an effective pore fluid, then Gassmann&#39;s fluid substitution is utilized to model the elastic properties. The effective pore fluid is modeled assuming Brie&#39;s empirical relationship (Brie et al. 1995). 
     Alternative Fluid Substitution Methodologies for 4D Studies 
     An alternative fourth approach to modeling the effective elastic properties of a fluid saturated rock is Gassmann-consistent and satisfies both the low and high frequencies limits on the elastic moduli. This allows pore fluids to be distributed throughout the rock in a realistic way, described by a correlation function. To model the elastic moduli, a spatial distribution of the pore fluid must be assumed and also a length scale which describes its spatial distribution. The method is described in Toms et al. (2006). This approach allows the elastic moduli to be modeled when immiscible pore fluids are spatially distributed in a complicated way. That is, pore fluid spatial distribution can be modeled assuming that fluids are spatially distributed according to a correlation function (with an associated correlation length). When approach No, 3 is assumed, wave attenuation can also be calculated. 
     There are classes of fluid substitution methods that allow additional physical effects of frequency-dependent P-wave attenuation and phase velocity to be modeled for partially saturated rock. These methods are consistent with the low and high frequency limits on the elastic moduli of partially saturated rock, and allow the effects of wave induced fluid flow to be accounted for. These methods do not appear to have been applied in any 4D time lapse studies. The current invention uses the frequency behavior of attenuation and velocity to derive the properties of spatial distribution of heterogeneities in pore fluid. 
     Attenuation and Velocity Dispersion 
     Two important characteristics of wave propagation are velocity and attenuation. To date, most geophysical technology utilized to image the subsurface relies fundamentally upon velocity, which is typically assumed to be independent of frequency (i.e. not dispersive). One exception is in the borehole environment where Stoneley wave attenuation is being utilized to estimate permeability. There are two primary reasons why attenuation is currently not exploited. First, a best practice for estimating or measuring attenuation does not exist, and there is still much debate on whether attenuation needs to be measured or whether it can be reliably extracted from the geophysical data via processing. The second reason is that attenuation data are difficult to interpret. Much of this difficulty stems from: (a) not understanding the physical mechanisms responsible for the loss of wave energy and (b) not having access to appropriate rock physics models. 
     The present invention addresses the second issue. That is, given a set of reliable attenuation and/or velocity dispersion measurements, how could or should they be interpreted? The issue of how attenuation is to be measured or extracted from geophysical data is a separate concern. Nevertheless, there are presently techniques to do this, and it may be expected that these methods will be improved in the future, and that developments in interpretation strategies will act as further motivation. Specifically, the invention interprets attenuation in terms of three parameters significant for characterization of subsurface fluid distribution. They are the fluid spatial distribution, the length scale and the relative percentage of fluids. Hence, implicit within the invention is the assumption that wave attenuation is caused by one physical mechanism due to the presence of two or more fluids. 
     In general, wave attenuation may be caused by many different physical mechanisms resulting from different characteristics of the rock and/or fluids. That is, wave attenuation may result from the presence of fractures within porous fluid saturated rock (Galvin et al. 2007). It may be the result of heterogeneities in rock porosity and permeability (Pride et al. 2004). It may arise due to the presence of small cracks on the pore scale interacting with pore fluids (Gist 1994). In short, there are many different physical mechanisms which have the potential to cause significant attenuation if the physical context arises. As such, attenuation is diagnostic of specific rock and/or fluid features. All of these mechanisms result in what is called intrinsic attenuation of wave energy. There are of course other reasons why waves suffer attenuation, such as geometrical spreading, mode conversion and reflection (scattering). The latter effects lead to what is termed apparent attenuation, in which there is no energy loss. The invention assumes that apparent attenuation has been separated from intrinsic attenuation. That is, the invention operates on intrinsic attenuation values to estimate characteristics of subsurface fluid distribution. 
     The Physical Mechanism 
     When compressional waves encounter two or more fluids within a porous rock, there can be significant attenuation and dispersion of the propagating waveform. This occurs because pressure gradients are induced between regions of the rock saturated by different fluids. The pressure gradients cause fluid flow relative to the rock frame, which leads to the viscous dissipation of wave energy. Thus, work is being done by the compressional wave on the fluids within the porous rock, and consequently the energy of the acoustic wave is correspondingly reduced. In simple analogous terms, one could imagine compressing a sponge partially saturated by water. The compression of the sponge causes water to spread to other regions of the sponge previously unsaturated. It is this kind of fluid movement relative to the framework of the sponge or rock, which leads to the loss of wave energy. 
     As a compressional wave cycles between compression and rarefaction, pressure gradients are induced in opposing directions through out the rock in accordance with the wave frequency. Fluid pressure equilibration will occur only for low enough wave frequencies when the characteristic length scale of the fluid heterogeneities is smaller than the fluid pressure diffusion length (Norris, 1993). The diffusion length scale is defined by 
                 λ   d     =         κ   ⁢           ⁢   N     ωη         ,         
where κ is rock permeability, η is fluid viscosity, ω is the angular wave frequency and N is a combination of poroelastic moduli. When fluid pressure equilibrates, it is analogous to having a homogeneous saturation of fluid. On the other hand, when wave frequencies are sufficiently high (or the heterogeneity length scale is greater than the diffusion length, there is not sufficient time for fluid pressure equilibration to occur. Hence, there are individual patches through out the rock having constant but different fluid pressures. It is at intermediate wave frequencies or heterogeneity length scales where wave attenuation and velocity dispersion may be significant.
 
     The practical implication of this physical mechanism is that compressional waves of any frequency may be affected by attenuation and velocity dispersion. That is, it is not necessarily the case that fluid pressure will be equilibrated during seismic frequencies as the heterogeneity length scale needs to be smaller than the diffusion length. One clear example is the influence of shallow gas (or gas chimneys) on seismic surveys, where significant losses in compressional wave energy have been observed in geophysical data. 
     Experimental Observation of Saturation Heterogeneities and Associated Effects 
     Fundamental to understanding compressional wave propagation through partially fluid saturated rock is the need to consider fluid distribution, in particular its heterogeneity length scale. It is the heterogeneity length scale relative to the diffusion length of the rock which influences the peak frequency of attenuation. Cadoret et al. (1990) showed experimentally that saturation heterogeneities greater than the diffusion length may cause significant attenuation and velocity dispersion of ultrasonic waves. In a series of experiments, the fluid distribution of a carbonate rock undergoing changes in water saturation was observed utilizing X-ray tomographic imaging; acoustic measurements were taken concurrently. The images show that different fluid saturation processes result in different fluid distributions for the same average water saturation. The acoustic measurements indicate the velocity and attenuation at certain water saturations is influenced by the fluid distribution pattern. While the images acquired during the experiment indicate fluid spatial distribution, only qualitative measures were utilized to infer that the heterogeneity length scale was greater than the diffusion length. 
     Toms et al. (2009) analyzed a series of X-ray tomographic images of partially water and gas saturated limestone undergoing changes in saturation. Statistical methods were utilized to characterize changes in the spatial distribution of fluids as water saturation was altered. The statistical measures varied significantly for small incremental increases in gas saturation at larger water saturation. When gas saturation was large, the statistical measure characterizing pore fluid spatial correlation did not change significantly, but only the characteristic length scale altered. By applying a specific method of fluid substitution, a technique was presented to forward model wave attenuation and velocity dispersion for the true fluid distribution patterns as shown in X-ray tomographic images. Hence, this shows directly that the spatial distribution (and its characteristic length scale) of fluid heterogneities may impact wave attenuation and velocity dispersion as water saturation is altered. 
     Frequency Dependent Attenuation Asymptotes 
     From a mathematical standpoint, an asymptote is a linear function that stands in a specific relationship to another function, such that the distance between the asymptote and the function approaches zero as they tend to infinity. Frequency dependent asymptotic behavior of attenuation is speculated to be linked to the spatial distribution of rock and fluid heterogeneities. For example, as wave frequency approaches zero an oblique frequency asymptote describes the low-frequency dependent behavior of attenuation. Similarly as wave frequency approaches infinity another oblique frequency asymptote describes the high-frequency attenuation. The zero and infinite frequency asymptotes represent two different physical scenarios. They represent the limit when fluid pressure is relaxed (i.e. the same fluid pressures through out the rock) and when fluid pressure is unrelaxed (different fluid pressures through out the rock). The peak frequency of attenuation occurs at the intersection frequency of the asymptotes, this is the frequency at which attenuation reaches a maximum. The specific intersection frequency depends on rock permeability and heterogeneity length scale. Gurevich and Lopatnikov (1995) show that the obliqueness of the attenuation asymptotes depends on the type of spatial variation of heterogeneities. 
     This invention employs a model-based predicted frequency behavior such as an asymptote to determine subsurface fluid distribution. It is believed that the significance of attenuation asymptotes or the frequency behavior at high and low frequencies has not been widely known and appreciated, and that there are no published reports that utilize attenuation asymptotes to determine rock and/or fluid distribution characteristics. 
     The present invention aims at identifying three aspects of subsurface fluid saturation utilizing geophysical data. The first objective of the invention is to determine the fluid spatial distribution. That is, whether fluid spatial distribution is homogeneous or heterogeneous. For the case where pore fluid distribution is heterogeneous, the invention seeks to identify potential spatial patterns.  FIG. 1  shows two examples of potential hydrocarbon spatial distribution within an otherwise water saturated rock:  110  a random distribution of hydrocarbon patches (denoted by the symbol x) and  111  a periodic distribution of hydrocarbon patches. Alternatively (not shown), a periodic spatial distribution might be regularly distributed layers, and a random spatial distribution might be randomly distributed layers. In any case, the spatial distribution may or may not be described by, for example, a correlation function or a linear path function. 
     The second objective of the invention is to identify the size of sub-wavelength saturation heterogeneities.  120  shows a spatial distribution of hydrocarbon patches in an otherwise water saturated porous rock.  121  shows the characteristic length of the hydrocarbon patches.  122  shows that the characteristic length of the hydrocarbon heterogeneities is much less than the P wavelength  123  of the geophysical data. The third objective of the invention is to identify the relative percentage of saturating fluids, that is, how much water or hydrocarbon is present.  130  shows the relative saturation of two different fluids, given by volume fractions v 1 , v 2 .  131  shows that the porosity of the rock is fully fluid saturated by two different fluids of volume fractions v 1 , v 2 . 
     The invention employs physical relationships which relate spatial distribution, patch size and percent saturation to P wave attenuation and velocity dispersion. Hence, when P wave attenuation and dispersion is negligible, the invention cannot be employed to estimate subsurface fluid distribution parameters. 
     A flow chart showing basic steps in one embodiment of the invention is shown in  FIG. 2 . The invention operates on attenuation and/or velocity dispersion data  210  that have been either measured directly using instruments or have been estimated using processing techniques from geophysical data. The invention can operate on attenuation and velocity data obtained from any type of geophysical data, such as seismic, cross well tomography, VSP, sonic logs, laboratory measurements on core plugs or from utilizing any combination of the above listed data types. 
     The geophysical data types span a large frequency range up to hundreds of Hz for seismic to 1 kHz for logging to 1 MHz for laboratory measurements. Presumably, the processing methods for converting geophysical data into attenuation and velocity dispersion information will be different from data type to data type. The reason for this is that the data types have very different frequency bandwidth, different acquisition geometries, different wavefield characteristics (transmission versus reflection) and different sized volumetric regions impacting wave propagation. Common to all the data processing methods and measuring techniques applied to the geophysical data is the need to obtain attenuation and velocity data with respect to wave frequency. That is, it is necessary to have reliable attenuation/velocity values and frequency values. Furthermore, data processing methods that assume specific frequency dependent attenuation behavior should be general enough to estimate attenuation when it exhibits a range of frequency dependencies. That is, methods that assume Q is independent of frequency may provide attenuation-frequency estimates that are biased by a constant attenuation assumption. 
     A first embodiment of the invention requires, in addition to processed geophysical data, rock and fluid property information  200 . Specifically, the additional data required is dry rock elastic properties: bulk K d  and shear moduli μ d , grain bulk modulus K g  and density ρ g , rock porosity φ and permeability κ. The fluid properties required are bulk modulus K f , density ρ f  and viscosity η f . Correct fluid bulk moduli are important for obtaining good estimates of subsurface fluid saturation parameters. It is necessary to have accurate bulk modulus of each fluid, as the contrast in elastic compliance between fluids significantly affects the magnitude of attenuation. Hence, knowledge of reservoir pressure and temperature may be required to refine fluid bulk moduli estimates utilizing principles of thermodynamics. 
     At step  220 , an appropriate rock physics model is chosen. The rock physics model must be able to relate frequency dependent attenuation and velocity dispersion to heterogeneous features of the fluid distribution. That is, the model needs to consider:
         a spatial distribution of two or more fluids, like periodic  111  or random  110  patches, to be selected by the user according to, for example, his/her experience in the area of what fluid distributions are more likely, or an empirical determination from saturation experiments on core plugs or fluids, or previous interpretation of similar studies, or the user may examine the asymptotic frequency dependency of the attenuation, among other possibilities for guessing a likely spatial distribution;   a characteristic length scale  121  of the heterogeneous distribution;   a relative percentage  130  of fluids, which means the volume fraction of fluid types within an effective volume. (Sometimes the fluid types within an effective volume may be only patches of either oil, gas or water; Alternatively, within the effective volume there may be fluid types within patches that are composites of two or more fluids, such as a mixture of oil and water or oil and gas)       

     The rock physics model should also include all necessary physical dependencies that impact frequency dependent wave attenuation and velocity dispersion. This means the rock physics model should include all physical factors that result in
         translational shifting of the attenuation and velocity dispersion curves in frequency; (see the examples discussed below in connection with  FIGS. 3-7 )   magnitude variation of the attenuation and velocity dispersion (see examples discussed below in connection with  FIGS. 8-9 ).       

     There are several different rock physics models that satisfy the above requirements. Some appropriate rock physics models include: White et al. 1977, Johnson 2001, Mueller and Gurevich 2004, Toms et al. 2007, and Mueller et al. 2008. There are other rock physics models that may also be utilized within this invention. To assist with choosing an appropriate rock physics model, the attenuation and velocity dispersion data  210  could be analyzed if multiple frequency values were obtained, i.e., if attenuation and velocity have the form {Q 1 ,ω 1 ; Q 2 ,ω 2 ; Q 3 ,ω 3 ; . . . } and {v 1 ; ω 1 ; v 2 ,ω 2 ; v 3 , ω 3 ; . . . }, respectively. The reason is that most rock physics models (with the exception of Mueller and Gurevich 2004 and Toms et al. 2007), specifically assume a spatial distribution of the pore fluids. The spatial distribution greatly impacts the frequency dependent attenuation characteristics predicted by the model {{tilde over (Q)} 1 ,{tilde over (ω)} 1 ; {tilde over (Q)} 2 ,{tilde over (ω)} 2 ; {tilde over (Q)} 3 ,{tilde over (ω)} 3 ;}. Hence, the observed data values {Q 1 ,ω 1 ; Q 2 ,ω 2 ; Q 3 ,ω 3 ; . . . } could be used to discriminate between rock physics model predictions {{tilde over (Q)} 1 ,{tilde over (ω)} 1 ; {tilde over (Q)} 2 ,{tilde over (ω)} 2 ; {tilde over (Q)} 3 ,{tilde over (ω)} 3 ;}, so as to eliminate models that do not predict similar behavior. 
     At step  230 , the frequency dependent asymptotes of the rock physics model are derived or already derived (known) asymptotes are obtained for the rock physics model. These are relationships that define the limit of attenuation behavior over a certain frequency range. Usually, two asymptotes can be derived from the rock physics model by considering the limiting cases of model behavior when ω→0 and ω→∞. For certain rock physics models, which have characteristic length scales and associated frequencies ω i  or ω 1  and ω 2 , it may also be possible to derive other asymptotes ω→ω i  or ω 1 &lt;ω&lt;ω 2  for certain model parameters. However, the two unique asymptotes ω→0 and ω→∞ implicitly contain necessary information on subsurface fluid distribution parameters and can be used together (without other asymptotes) to interpret the attenuation behavior in  210 . In this invention, it is preferred to use at least two attenuation asymptotes for any given model. 
     For example, consider an embodiment of the invention where the rock physics model of Mueller and Gurevich (2004) is used. If the correlation function that describes the spatial distribution is given by periodically alternating layers, the low frequency asymptote of the attenuation is 
                       Q     -   1       =         S   1     ⁢     S   2     ⁢     L   H   2     ⁢   s   ⁢           ⁢   ω       12   ⁢     D   0           ,           (   3.1   )               
where S 1 , S 2  refer to the relative volume fraction of fluids, L H  is the heterogeneity length scale, s a function of the low and high frequency elastic moduli, D 0  is the diffusitivity and ω is the frequency. The high frequency asymptote for attenuation is
 
                       Q     -   1       =     1     1   +         S   1     ⁢     S   2     ⁢     L   H     ⁢       2   ⁢   ω           δ   ⁢       D   0                 ;           (   3.2   )               
here δ is a function of the low and high frequency elastic moduli. The asymptotes (3.1) and (3.2) derived for this model contain equivalent information on subsurface fluid distribution as the model of Müller and Gurevich (2004) with periodic correlation function.  FIG. 11  shows attenuation (Q −1 ) versus wave frequency predicted using the rock physics model ( 1150 ), and attenuation estimated using the low ( 1060 ) and high ( 1070 ) frequency asymptotes. Away from the frequency of peak attenuation ( 1155 ), the asymptotes and model are equivalent. In the vicinity of the frequency of peak attenuation, hereinafter also called “peak frequency,” the asymptotes do not predict the attenuation correctly; however, they can still be utilized in interpreting attenuation around the peak frequency (step  240 ). This is demonstrated in the example section. The reason is the two asymptotes are unique for each rock physics model, such that attenuation variation with frequency is uniquely defined once two asymptotes are known.
 
     The derived or known asymptotes of a specific rock physics model can assist with determining whether the rock physics model employed to interpret the attenuation is appropriate. That is, if the attenuation asymptotes do not exhibit the same frequency dependence as the geophysical data {Q 1 ,ω 1 ; Q 2 ,ω 2 ; Q 3 ,ω 3 ; . . . } then that specific rock physics model should not be chosen. The exception is when the attenuation data {Q 1 ,ω 1 ; Q 2 ,ω 2 ; Q 3 ,ω 3 ; . . . } is constant with frequency. Then it is not possible to use the asymptotes to evaluate whether a specific rock physics model should be chosen or not. Where there is a good correspondence between asymptotes and the geophysical data, it indicates the likely spatial distribution of heterogeneities, like periodic  110  or random  111  patches, or random and periodic layering (see  FIG. 10 ). 
     At step  240 , each asymptote is inverted, so that one subsurface parameter is a function of attenuation, frequency and the second subsurface parameter. In general, the attenuation asymptote will take the form
 
 Q   −1   =f (ω, L   H   Q   −1   ,S   i   ,P ),  (3.3)
 
where ω is frequency, L H  refers to the heterogeneity length scale ( 120 ), S i  refers to the relative volume fraction ( 130 ) of fluids and P is a parameter containing rock and fluid property information. This function can be inverted so as that the heterogeneity length scale L H  is a function of the other parameters:
 
 L   H   =f (ω, Q   −1   ,S   i   ,P ).  (3.4)
 
     For instance, assuming a periodic distribution of fluids, the low frequency attenuation asymptote (3.1) is given by 
                     L   H     =       [       12   ⁢     D   0         ω   ⁢           ⁢       S   1     ⁡     (     1   -     S   1       )       ⁢   sQ       ]       1   2               (   3.5   )               
and the high frequency attenuation asymptote (3.2) is
 
     
       
         
           
             
               
                 
                   
                     L 
                     H 
                   
                   = 
                   
                     
                       
                         
                           D 
                           0 
                         
                         
                           2 
                           ⁢ 
                           ω 
                         
                       
                     
                     ⁢ 
                     
                       
                         
                           
                             ( 
                             
                               Q 
                               - 
                               1 
                             
                             ) 
                           
                           ⁢ 
                           δ 
                         
                         
                           
                             S 
                             1 
                           
                           ⁡ 
                           
                             ( 
                             
                               1 
                               - 
                               
                                 S 
                                 1 
                               
                             
                             ) 
                           
                         
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   3.6 
                   ) 
                 
               
             
           
         
       
     
     Other fluid distributions would have different forms of asymptotes (3.3) which would result in different length scale-attenuation relationships of which (3.4) is the generalized functional form. It is preferred to invert at least two frequency asymptotes for any given rock physics model, that is, a low and a high frequency asymptote. In Table 1, asymptotes derived from specific rock physics models are given. These are some of the inverted asymptotes that may also be utilized to interpret the geophysical data. In Table 1, S 1 , S 2  refer to the relative volume fraction of fluids, L H  is the heterogeneity length scale, s a function of the low and high frequency elastic moduli, D 0  is the diffusivity coefficient and ω is the frequency. Alternatively, the relative volume fraction v i  may be expressed as a function of the other parameters such that
 
 S   i   =f (ω, Q   −1   ,L   H   ,P *).  (3.7)
 
     However, using the latter approach is difficult as the P parameters also have dependencies on the relative volume fractions v i  which need to be isolated to form P*. 
     
       
         
           
               
             
               
                 TABLE 1 
               
             
            
               
                   
               
               
                 Displaying inversion scenarios for the saturation dependencies 
               
            
           
           
               
               
               
               
            
               
                   
                   
                   
                 Saturation 
               
               
                   
                   
                   
                 heterogeneity length 
               
               
                   
                   
                 Q −1  frequency 
                 scale L H  from 
               
               
                   
                 Geometry 
                 dependence 
                 asymptote 
               
               
                   
                   
               
               
                   
                 Periodic layering low-frequency asymptote 
                 ω 
                 
                   
                     
                       
                         
                           [ 
                           
                             
                               12 
                               ⁢ 
                               
                                 D 
                                 0 
                               
                             
                             
                               ω 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 
                                   S 
                                   1 
                                 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     1 
                                     - 
                                     
                                       S 
                                       1 
                                     
                                   
                                   ) 
                                 
                               
                               ⁢ 
                               sQ 
                             
                           
                           ] 
                         
                         
                           1 
                           2 
                         
                       
                     
                   
                 
               
               
                   
                   
               
               
                   
                 Random layering low-frequency asymptote 
                 {square root over (ω)} 
                 
                   
                     
                       
                         
                           
                             
                               2 
                               ⁢ 
                               
                                 D 
                                 0 
                               
                             
                             ω 
                           
                         
                         ⁢ 
                         
                           1 
                           
                             s 
                             ⁡ 
                             
                               ( 
                               
                                 Q 
                                 + 
                                 1 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                 
               
               
                   
                   
               
               
                   
                 Periodic layering high-frequency asymptote 
                 
                   
                     
                       
                         1 
                         
                           ω 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             
                               D 
                               0 
                             
                             
                               2 
                               ⁢ 
                               ω 
                             
                           
                         
                         ⁢ 
                         
                           
                             
                               ( 
                               
                                 Q 
                                 - 
                                 1 
                               
                               ) 
                             
                             ⁢ 
                             δ 
                           
                           
                             
                               S 
                               1 
                             
                             ⁡ 
                             
                               ( 
                               
                                 1 
                                 - 
                                 
                                   S 
                                   1 
                                 
                               
                               ) 
                             
                           
                         
                       
                     
                   
                 
               
               
                   
                   
               
               
                   
                 Random layering high-frequency asymptote 
                 
                   
                     
                       
                         1 
                         
                           ω 
                         
                       
                     
                   
                 
                 
                   
                     
                       
                         
                           
                             
                               D 
                               0 
                             
                             
                               2 
                               ⁢ 
                               ω 
                             
                           
                         
                         ⁢ 
                         
                           ( 
                           
                             1 
                             - 
                             Q 
                           
                           ) 
                         
                         ⁢ 
                         δ 
                       
                     
                   
                 
               
               
                   
                   
               
            
           
         
       
     
     At step  250 , the inverted asymptotes are utilized to interpret geophysical attenuation data {Q 1 ,ω 1 ; Q 2 ,ω 2 ; Q 3 ,ω 3 ; . . . }; velocity dispersion data can also be utilized in this step {v 1 ; ω 1 ; v 2 ,ω 2 ; v 3 , ω 3 ; . . . }. The approach is to substitute one attenuation value {Q 1 ,ω 1 } into each inverted asymptote. Then it is necessary to iterate through the parameter space of the independent subsurface parameter in a physically reasonable way. That is, if heterogeneity length scale L H  is functionally dependent on relative fluid saturation v i  (i.e. inverted asymptote (3.4)), it is necessary to vary relative fluid saturation from 0% to 100% (see  FIGS. 13-16 ). This must be conducted for each inverted asymptote. The intersection of the inverted asymptotes uniquely determines the heterogeneity length scale L H  ( 270 ) and the relative fluid saturation v i  ( 280 ) of most subsurface fluid distributions. To obtain the spatial distribution of different pore fluids ( 260 ), one may try different fluid distributions to estimate the length scales and saturation. The fluid distribution (i.e. the spatial distribution of different pore fluids) constrains the frequency-dependent attenuation behavior. Most rock physics models are not general and are derived for specific spatial distributions of pore fluids. In these cases, the attenuation and velocity relationships are valid only for the assumed fluid distribution. However, there are several rock physics models (e.g., Gurevich and Lopatinkov 1995, Müller and Gurevich 2004, Toms et al 2007) that are general in that the attenuation and velocity frequency behavior can be specified in terms of a function of the spatial distribution of porefluids, and hence spatial distribution is a variable, i.e a parameter, of the model. However, in practice one then assumes a distribution in terms of a correlation function to estimate what the attenuation/velocity behavior would be. If there are multiple attenuation-frequency points it may be possible to eliminate certain fluid distributions at the outset because of the asymptotic behavior. 
     There may be benefits to using a general rock physics model for the analysis, even though one has to choose a distribution. Within generalized rock physics models, the mathematical formulation/framework for each different distribution is the same. When one switches between different rock physics models, which a priori assume specific fluid distributions, the underlying mathematical formulation for the poro-elastic parameters can be very different. This complicates the use of inversion to extract saturation and length scale information. 
     On occasion there may be two intersection values of the heterogeneity length scale L H  and relative fluid saturation v i . In these cases, velocity dispersion could be used as extra information to determine the correct fluid saturation. Alternatively, the values can be re-substituted into the full rock physics model as described in Example 4 and the best results chosen based on the behavior of the rock physics relationships. In circumstances where two intersection points occur, the full rock physics model needs to be employed utilizing the heterogeneity length scale and relative fluid saturation relationships and {Q 1 ,ω 1 ; Q 2 ,ω 2 ; Q 3 ,ω 3 ; . . . } and velocity {v 1 ; ω 1 ; v 2 ,ω 2 ; v 3 , ω 3 ; . . . } data. By comparison of either geophysical data set with the model responses utilizing the heterogeneity length scale and saturation relations a unique intersection point can be interpreted, which indicates the actual length scale and relative saturation. 
     The description above is the forward-modeled utilization of the inverted asymptotes (which works with attenuation and velocity values at one or more frequencies) to determine the subsurface fluid distribution parameters. With multiple measurements of attenuation and velocity as a function of frequency, the inversion can be set up using a forward model of attenuation and dispersion based on the rock physics relationships and asymptotes and solving for the saturation parameters with the greatest likelihood of matching the data in a least squares sense. The analysis may be done first for individual frequency measurements (see  FIG. 13 ). Then an average of these results may be used as initial model and the attenuation and dispersion predicted using the rock physics relationship (see  FIGS. 14-16 ). The model may then be updated so that the predicted values better match the measurements within the constraints of reasonable physical behavior and geological knowledge. Preferably, this inversion scheme would weight the attenuation and dispersion values by their relative uncertainties allowing for the fact that measurements at one frequency range (for example, for sonic logging frequencies) might be more accurate than measurements at another frequency range (for example, for seismic frequencies). 
     EXAMPLES 
     Compressional wave attenuation and velocity dispersion are affected by many different characteristics of pore fluid saturation and rock properties. In the examples shown below three different behaviors are examined. The first example demonstrates the influence of rock permeability and heterogeneity length scale on wave attenuation and velocity dispersion. It is shown that changes in rock permeability or heterogeneity length scale result in lateral shifts of attenuation/dispersion curves with frequency. The second set of examples explores the impact that rock stiffness and fluid contrast have on the magnitude of attenuation. The influence of average fluid saturation on attenuation is a secondary factor that also affects the magnitude of attenuation. The third example illustrates the influence of different fluid spatial distributions on wave attenuation and velocity dispersion. It is demonstrated that random and periodic layering of subsurface fluid distributions results in different frequency dependency of the attenuation. The last set of examples illustrates a workflow for prediction of subsurface fluid distributions from attenuation. 
     Example 1 
     Where is the Peak Frequency of Attenuation? 
     In this example, P-wave velocity dispersion and attenuation are modeled assuming a periodic distribution of light oil and water with a porosity 22%. The Biot coefficient of the rock is 0.824 (see definition below). Two different cases are examined. In  FIGS. 3A-3C and 4A-4C , the characteristic length scale of the spatial distribution of fluid is fixed at 40 cm, while rock permeability is altered from 1 Darcy to 10 milli-Darcy. In  FIGS. 5A-5C and 6A-6C , the characteristic length scale of the fluid distribution is altered from 4 m to 4 cm, while rock permeability is fixed at 1 Darcy. The attenuation curves (e.g.,  430  in  FIGS. 4A-4C ) are plotted in logarithmic inverse Q to demonstrate the frequency dependence of attenuation. Other representations of Q with respect to log frequency, such as Q or inverse Q do not reveal the exact frequency dependence of the attenuation. Similarly, the velocity dispersion curves (e.g.  330  in  FIGS. 3A-3C ) are plotted on a logarithmic frequency scale. In each drawing, the shaded regions ranging from lightest to darkest correspond to different geophysical observation bandwidths: seismic (e.g.,  340 ), VSP ( 350 ), cross-well ( 360 ), borehole logging ( 370 ), and laboratory ( 380 ), respectively. 
     When rock permeability changes from 1 Darcy to 10 milli-Darcy ( FIGS. 3A-3C ), the dispersion curve for velocity is shifted to lower frequencies. The translation across frequency influences the velocity dispersion observed within each bandwidth. For instance, across the seismic bandwidth [10-60 Hz] there is approximately 1% (˜25 m/s), 3% (100 m/s) and 1% (25 m/s) velocity dispersion estimated for each scenario. In addition, the functional dependence of velocity on frequency varies from being nonlinear to linear to nonlinear. Within the seismic bandwidth, the most permeable rock ( 3 A) has the lowest P-wave velocities, whilst the least permeable rock ( 3 C) has the highest wave velocities. Similar observations can be made for the other geophysical bandwidths. 
       FIGS. 4A-4C  show that changes in rock permeability from 1 Darcy to 10 milli-Darcy also result in lateral shifts of the attenuation curve  430  with respect to frequency. This can easily be observed by considering the point at which attenuation is a maximum, the corresponding frequency is known as the peak frequency of attenuation. The peak frequency of attenuation also represents the frequency at which P wave velocity changes most rapidly with respect to frequency. As rock permeability decreases the peak frequency of attenuation occurs within different bandwidths from cross well to seismic to low frequency (1 Hz). The translation of the attenuation curve with respect to frequency has a two fold impact. Firstly the magnitude of attenuation within each bandwidth changes and secondly the frequency dependent behavior changes significantly. For instance, within the seismic bandwidth from 10 Hz to 60 Hz, the attenuation behavior increases linearly with frequency, is constant with frequency or decreases with square root of frequency. The translation also impacts the magnitude of attenuation at each frequency, the greatest magnitude of attenuation with the seismic bandwidth occurs when the peak frequency is centered within this bandwidth. When the peak frequency is centered at other bandwidths, the magnitude of attenuation is decreased; although it still is significant if the peak frequency is situated within a neighboring frequency bandwidth. 
     It may be noted from  FIGS. 3A-3C and 4A-4C  that peak frequency is strongly influenced by the permeability. Therefore, the present inventive method should also be suitable for estimating permeability as well as, among others, fluid saturation and heterogeneity length scale. 
     In  FIGS. 5A-5C and 6A-6C , rock permeability is fixed at 1 Darcy, but the characteristic length scale of the light oil patches is altered.  FIGS. 5A-5C  show that altering the characteristic length scale of the heterogeneous patches of oil causes the velocity dispersion curve to translate with respect to frequency. That is, as the characteristic length scale varies from 4 m to 4 cm the dispersion curve is shifted to higher frequencies. Translation of the dispersion curve with respect to frequency alters both the measured P wave velocity within each bandwidth, the gross percentage of dispersion and the velocity dependence on frequency. For instance, within the seismic bandwidth, when the patch size is 4 m the velocity measured varies from 3500 m/s to 3525 m/s, which represents less than 1% velocity dispersion. When the heterogeneity length scale is 40 cm the velocity decreases rapidly to 3325-3300 m/s, with less than 1% velocity dispersion. Finally, when the characteristic length scale is 4 cm, the velocity is constant within the seismic bandwidth at 3300 m/s. 
     Similarly,  FIGS. 6A-6C  show that altering the characteristic length scale of the light oil heterogeneities also causes lateral shifting of the attenuation curve with respect to frequency. As the characteristic length scale of the heterogeneities decreases from 4 m to 4 cm, the peak frequency of attenuation shifts to higher frequencies. Frequency translations of the attenuation curve result in changes in attenuation magnitude and frequency dependent behavior within a given bandwidth. In this scenario, for the seismic bandwidth (10-60 Hz), 4 m heterogeneities cause attenuation to decrease with frequency, while 40 cm and 4 cm heterogeneities cause attenuation to increase with frequency. 
     When variations in rock permeability are considered in conjunction with alterations in the characteristic length scale of fluid heterogeneities, it is possible to determine the peak frequency location.  FIG. 7  shows a crossplot of permeability and characteristic length scale. Each quantity is plotted in log base 10. This means that for the heterogeneity length scale, 0 corresponds to 1 m, −1 corresponds to 10 cm, −2 corresponds to 1 cm, −3 corresponds to 1 mm, etc. The gray-scale bar represents the peak frequency of attenuation. Light colors indicate that the peak frequency occurs at higher frequencies, while darker colors indicate that the peak frequency occurs at low frequencies. Three different lines are shown in the diagram according to different geophysical measurements.  730  shows that there is a range of rock permeabilities and characteristic length scales for which attenuation peaks within the seismic bandwidth. Similarly, line  720  shows attenuation will peak within the logging frequencies for a different set of rock permeabilities and characteristic length scales, while line  710  indicates that the peak frequency will be located at even higher frequencies for the laboratory frequency band. 
     To summarize, alterations of rock permeability or characteristic length scale result in lateral frequency shifts of attenuation and dispersion. Alterations in rock permeability or characteristic length scale do not change the peak magnitude of attenuation (at the peak frequency) or the gross percentage of velocity dispersion. However, when attenuation and dispersion are analyzed within a single observation bandwidth, the resulting lateral shift of attenuation and dispersion curves result in changes of attenuation magnitude/behavior and percentage dispersion across the observation bandwidth. When permeability and length scale changes are considered together, it is possible to define the location of the peak frequency of attenuation in terms for different permeable rocks and different heterogeneity length scales. It is evident that there is a large set of possible permeability rocks, which have characteristic heterogeneity length scale which result in the peak frequency of attenuation being situated in any one observation bandwidth. 
     Example 2 
     What is the Magnitude of Attenuation at the Peak Frequency? 
     In this example the magnitude of attenuation is analyzed as a function of rock frame stiffness and fluid bulk moduli contrast. Both factors significantly impact the magnitude of attenuation across the entire frequency bandwidth. Specifically, the magnitude of attenuation at the peak frequency will be examined. This represents the maximum attenuation that can arise for a specific combination (and concentration) of pore fluids within a rock of certain elastic stiffness. Here the rock frame stiffness is parameterized in terms of the Biot coefficient α R , which is a measure of bulk macroscopic compressibility of a volume of dry rock relative to the compressibility of its constituent materials. The Biot coefficient is defined as 
                 α   R     =     1   -       K   d       K   g           ,         
where K d  is the bulk modulus of the dry rock and K g  is the bulk modulus of the grain material. When α R →0.5, the rock is stiff and the acoustic response (attenuation, velocity) of a partially saturated rock will be dominated by the acoustic characteristics of the dry rock. When α R →1.0, the rock is very compressible and the acoustic response (attenuation, velocity) will be significantly impacted by the saturating fluid. Both the relative concentration of fluids and the contrast in fluid bulk moduli between respective pore fluids may have significant impact on the acoustic properties (attenuation, velocity). Here the fluid contrast α F  is defined as
 
                 α   F     =     1   -       K     f   ⁢           ⁢   1         K     f   ⁢           ⁢   2             ,         
where K f1  is the bulk modulus of the most compressible fluid and K f2  is the bulk modulus of the stiff (least compressible) fluid. When α F →1, fluid bulk moduli differ significantly, hence there is large contrast in fluid compressibilities. However, when α F →0.5, the compressibility of the pore fluids are of the same order of magnitude and there is little contrast in their elastic properties. The elastic properties of fluids such as water, oil and gas depend significantly on subsurface pressure and temperature. Also, the bulk moduli of water, oil and gas is affected by the specific chemical compositions, respectively, i.e. salt, API, etc.
 
     In  FIG. 8 , the Biot coefficient of the rock is altered from α R =0.5→1.0 and the fluid contrast is altered from α F =0.5→1.0. The peak magnitude of attenuation given by inverse Q is shown in the color bar as log base 10. The darkest colors correspond to small attenuation, while the lightest colors indicate significant attenuation. The data values above line  1010  indicate regions where attenuation is significant. The relative concentration of the most compressible fluid to the least compressible fluid is v f1 =0.10. The spatial distribution of the pore fluids is described by an exponential correlation function. The fluid bulk modulus of the least compressible fluid is taken to be K f2 =2.84 GPa corresponding to water. The grain bulk modulus is K g =76 GPa, which corresponds to the elastic moduli of calcite grains. The horizontal line  1020  indicates likely gas-water response, while  1030  indicates likely light oil-water response and line  1040  indicates higher API oil and water. For any given line, the attenuation decreases as rock stiffness increases (i.e. traversing right to left across the picture). Clearly, the greater the contrast in fluid bulk moduli (line  1020  versus line  1040 ), the greater the potential attenuation for any given rock stiffness. 
     In  FIG. 9 , the relative concentration of the most compressible fluid to the least compressible fluid is v f1 =0.40. There is a slight change in the magnitude of attenuation (shades of gray) across the parameter space. In general there is increase in attenuation for the lower contrast fluid systems ( 1130 - 1140 ), while a reduction in attenuation observed for the largest contrast systems ( 1120 ). The data values above line  1110  indicate regions where attenuation is significant. Line  1120  represents attenuation values for a range of rock compressibilities for a system of 40% gas with 60% water. Line  1130  shows attenuation values for a range of rock compressibilities with 40% light oil and 60% water. As the concentration of less compressible fluids increases attenuation is significant over a larger range of rock compressibilities. For instance, line  1140  indicates that attenuation is significant for certain rocks even for a heavier oil water system. 
     Example 3 
     The Influence of Spatial Distribution on Frequency Dependent Attenuation 
       FIGS. 10A-B  demonstrate the effect of different spatial distributions on the frequency dependent characteristics of attenuation ( 10 A) and velocity dispersion ( 10 B). While only random and periodic layering of fluid heterogeneities were modeled, the effect of having different spatial distributions is well documented in the literature; see for example Mueller et al. 2008. The impact is that attenuation will increase or decrease more quickly for some fluid distribution patterns than others. 
     In  FIGS. 10  A-B, attenuation and velocity dispersion is estimated using the Mueller and Gurevich (2004) rock physics model by assuming two different spatial distributions of fluids. Here  1050  and  1051  refer to attenuation and velocity dispersion estimated assuming a random distribution of heterogeneous layers. Here  1060  and  1061  refer to attenuation and velocity dispersion estimated using a periodic distribution of heterogeneous layers. The difference in how attenuation and velocity varies with frequency as frequency decreases, arises because of the different spatial distributions of pore fluids. If measured geophysical values are present at different frequencies, like attenuation ( 1070 ,  1075 ,  1080 ) and velocity ( 1071 ,  1076 ,  1081 ) then the frequency dependence of the data points can be used to assist with determining the spatial distribution of pore fluids (and hence appropriate rock physics models to use in the invention). In this case, it is obvious that random layering of fluid heterogeneities produces frequency dependent behavior which is different from the geophysical data points, however models which assume periodic layering could potentially be utilized to interpret the data.  1041  and  1031  represent the velocity theoretical lower and upper limits of velocity which are independent of spatial distribution of fluids. 
     Example 4 
     Figures Demonstrating the Method for Interpretation of Intrinsic Q 
       FIG. 12  demonstrates the key features of utilizing the asymptotes to interpret attenuation behavior. The attenuation  1210  is modeled using the rock physics model of Mueller and Gurevich 2004 assuming periodically alternating layers. Both the low ( 1250 ) and high ( 1260 ) frequency asymptotes are given in equations (3.1) and (3.2). Attenuation data points are selected across the spectrum of frequencies (low:  1220 , peak:  1240  and high:  1230 ). In reality these data points would be supplied from the geophysical data. In  FIG. 13 , the inverted asymptotes (Eqs. (3.5) and (3.6)) are plotted as a function of water saturation utilizing one attenuation and frequency point.  1320  is the inversion of the low-frequency asymptote  1250 , and  1330  is the inversion of the high-frequency asymptote  1260 . Referring to  FIG. 13 , the intersection ( 1340 ) of the two inverted asymptotes ( 1320  and  1330 ) defines the length scale and relative saturation of the fluid distribution. Irrespective of whether observation point  1220 ,  1230  or  1240  is substituted into the inverted asymptotes, all produce a consistent length scale-saturation relation; i.e. all produce the same low ( 1320 ) and high ( 1330 ) frequency inverted asymptotes. Thus, interpretation can be performed with a single attenuation or frequency pair, regardless of lack of knowledge about whether it is on the low or high frequency side of the peak. With measurements at multiple frequencies, we can reduce uncertainties associated with measurement error by separately analyze each pair and averaging the results or solving for the most likely parameters using least squares. In addition to using attenuation, if available, velocity dispersion can also be analyzed to independently determine water saturation. 
     On the other hand, an issue with using, for example, seismic data and also log data, aside from REV (representative elementary volume) fluid distribution concern, is that often the attenuation values are only significant within a certain window of observation, as several of the figures show, and hence from that standpoint it does not make sense to combine data where measurement errors may be considerably different in an attempt to add additional certainty to the inversion. Moreover, one does not get additional physical information from adding more frequency points, but rather from utilizing different asymptotes. Whether data redundancy helps or not is a decision the user of the invention can make in the context of the particular application. 
     Further analysis using the rock physics relationships may also be useful. In  FIG. 14-16 , substitution of the length scale-saturation relations back into the rock physics model allows one to determine whether the frequency point  1220  (or  1230  or  1240 ) is located at the peak frequency of attenuation or at frequencies above or below the peak frequency.  1420  is the rock physics model constrained by the inverted asymptote  1320 .  1430  is the rock physics model constrained by the inverted asymptote  1330 . ( 1440  is their intersection, which occurs at the same value of water saturation as intersection  1340  in  FIG. 13 .) For example, curve  1330  may be obtained by substituting the frequency/attenuation data point  1220  into an inverted high-frequency asymptote equation such as equation (3.6), and then using that relationship to plot heterogeneity length scale versus saturation. Experience teaches that the constrained curve that is flattest with respect to variation in attenuation as water saturation changes depicts the frequency location of the attenuation point. That is, in this example the frequency point ( 1220 ) is located on the low frequency asymptote as the constrained rock physic model ( 1420 ) is flat.  FIG. 15  shows the same constrained rock physics curve using the input data from  1230 ; however, this time the input data frequency point is shown to be located on the high frequency asymptote  1260 , as the curve  1560  is flattest. ( 1550  is the inverted low-frequency asymptote and  1570  is the intersection point.)  FIG. 16 , which uses the frequency point  1240 , shows both constrained curves  1680  and  1690  vary significantly with water saturation. This indicates that the frequency point  1240  occurs in the vicinity of the peak frequency. The intersection points  1440 ,  1570 , and  1695  indicate relative fluid saturation. 
     Example 4 shows that one can determine where the Q value at a specific frequency is relative to the peak frequency of attenuation. That is, one can work out whether its frequency is less than the peak frequency and hence rests on the low frequency asymptote or occurs at the peak frequency of attenuation (in which case it does not rest on either asymptote) or whether the value rests on the high frequency asymptote. The benefit of knowing this information is that one can determine whether attenuation is increasing, decreasing, or is approximately constant with frequency. The above-described re-substitution technique of example 4 is especially useful when the curves in  FIG. 13  cross in two places. In that event, the asymptotes will change rapidly in the vicinity of saturation for one of the intersection points, which is not correct, but will be constant in the vicinity of the saturation, which is correct, and hence one can decide between the two points. 
     It is believed that the advantages of the invention have been explained in the foregoing disclosure, but to summarize a key point, the invention is important because otherwise saturation can never be estimated from attenuation in a physically meaningful way. That is, the interpretation will be biased by the choice of heterogeneity length scale. The invention estimates this length scale in addition to the saturation. Given multiple frequency data points, it is also possible to define the likely spatial distribution of the pore fluids, e.g. whether it is stratified or fractal, using the invention. The previously described patent publication by Hawthorne et al. primarily seeks fluid type rather than average saturation, and does not seek to identify length scale or spatial distributions of fluid. Tserkovank and Johnson (2001) do not employ asymptotic behavior of a rock physics model. Unlike the present invention, they do not solve for fluid saturation along with other fluid distribution parameters (spatial distribution, length scale and average saturation). 
     The foregoing patent application is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. For example, the invention is described mostly in terms of seismic attenuation and geophysical acquisition geometry of a land or marine survey. Seismic waves have a limited frequency bandwidth generally about 13-30 Hz. The method could alternatively be applied to any compressional wave (P-wave) attenuation measured or extracted from, for example, borehole, vertical seismic profile (“VSP”), or cross-well tomography, i.e. different frequency bandwidths and acquisition geometries. All such modifications and variations are intended to be within the scope of the present invention, as defined in the appended claims. Persons skilled in the art will readily recognize that in practical applications of the invention, at least some of the steps in the present inventive method are performed on or with the aid of a computer, i.e. the invention is computer implemented. 
     REFERENCES 
     
         
         Brie, A., Pampuri, F., Marsala, A. F. and Meazza, O., “Shear sonic interpretation in gas bearing sands,” SPE 30595 , Society of Petroleum Engineers Annual Technical Conference and Exhibition, Dallas , (1995). 
         Cadoret, T., Marion, D. and Zinszner, B., “Influence of frequency and fluid distribution on elastic wave velocities in partially saturated limestones,”  J. of Geophysical Research  100, 9789-9803 (1995). 
         Cadoret, Mavko, and Zinszner, “I Fluid distribution effect on sonic attenuation in partially saturated limestones,”  Geophysics  63, 154-160 (1998). 
         Galvin, R. J., Toms, J. and Gurevich, B., “A simple approximation for seismic attenuation and dispersion in a fluid-saturated porous rock with aligned fractures,”  SEG Expanded Abstracts  26, 1679 (2007). 
         Gassmann, F., “Elastic waves through a packing of spheres,”  Geophysics  16, 673-685 (1951). 
         Gist, G. A., “Interpreting laboratory velocity measurements in partially gas-saturated rocks,”  Geophysics  59, 1100-1109 (1994).” 
         Gurevich, B. and Lopatnikov, S. L., “Velocity and attenuation of elastic waves in finely layered porous rocks,”  Geophysical J. Int&#39;l.  121, 933-947 (1995). 
         Hill, R., “Elastic properties of reinforced solids: some theoretical principles,”  J. of the Mechanics and Physics of Solids  11, 357-372 (1963). 
         Johnson, D. L., “Theory of frequency-dependent acoustics in patchy-saturated porous media,”  J.l of the Acoustical Society of America  110, 682-694 (2001). 
         Müller, T. M. and Gurevich, B., “One-dimensional random patchy saturation model for velocity and attenuation in porous rocks,”  Geophysics  69, pp. 1166-1172 (2004). 
         Müller, T. M., Toms-Stewart, J., and Wenzlau, F., “Velocity-saturation relation for partially saturated rocks with fractal pore fluid distribution,”  Geophysical Research Letters  35, L09306 (2008). 
         Pride, S. R., Berryman, J. G. and Harris, J. M., “Seismic attenuation due to wave-induced flow,”  J. of Geophysical Research  109, No. B1, B01201 (2004). 
         Toms, J., Müller, T. M., Ciz, R. and Gurevich, B., “Comparative review of theoretical models for elastic wave attenuation and dispersion in partially saturated rocks,”  Soil Dynamics and Earthquake Engineering  26, 548-565 (2006). 
         Toms, J., Müller, T. M. and Gurevich, B., “Seismic attenuation in porous rocks with random patchy saturation,”  Geophysical Prospecting  55, 671-678 (2007). 
         Toms-Stewart, J., Müller, T. M., Gurevich, B., and Paterson, L., “Statistical characterization of gas patch distributions in partially saturated rocks,”  Geophysics  74, WA51-WA64 (2009). 
         Tserkovnyak, Y. and Johnson, D. L, “Can one hear the shape of a saturation patch,”  Geophysical Research Letters  29, 1108-1112 (2002). 
         White, J. E., “Computed seismic speeds and attenuation in rocks with partial gas saturation,”  Geophysics  40, 224-232 (1975). 
         White, J. E., Mikhaylova, N. G. and Lyakhovitskiy, F. M., ‘Low-frequency seismic waves in fluid-saturated layered rocks,”  Physics of the Solid Earth  (Transl. Izv.) 11, 654-659 (1976).