Inversion of 4d seismic data

CROSS-REFERENCE TO RELATED APPLICATIONS

0001This application claims the benefit of <patcit id="pcit0001" dnum="US92848007P" dnum-type="L"><text>U.S. Provisional application 60/928,480 which was filed on May 9, 2007</text></patcit>.

FIELD OF THE INVENTION

0002This invention relates generally to the field of hydrocarbon reservoir analysis, and more particularly to time-lapse seismic analysis. Specifically, the invention is a method for inversion of 4D seismic data to predict fluid and pressure change.

BACKGROUND OF THE INVENTION

00034D seismic, sometimes called time-lapse seismic, is a term that refers to performing an initial, or <i>base,</i> seismic survey of a subsurface region, then performing at least one later, or <i>monitor</i>, survey of the same region, attempting to duplicate the acquisition parameters and conditions as closely as possible. This enables comparison of seismic images of the subsurface over intervals of months or years of time to detect changes due to production of hydrocarbons or injection of fluids. More generally, the term <i>base</i> survey may refer to any survey performed earlier in time than the survey referred to as the <i>monitor</i> survey; i.e., the base survey may refer to an earlier monitor survey. 4D seismic is the only field-wide history match constraint for theoretical predictions of such changes by elaborate software modeling programs called reservoir simulators. Currently, one compares simulator results and 4D seismic data by comparing synthetic seismic data modeled using the simulator to measured 4D seismic difference data; alternatively, acoustic impedance (the product of density and acoustic velocity) modeled using a simulator is compared to inverted impedance, i.e. impedance (as a function of subsurface spatial location) inferred by <i>inversion</i> of seismic data, if inversion has been done. Ideally, the simulator model will be adjusted, or <i>updated</i>, using reservoir properties derived from quantitative 4D, i.e. saturation change and pressure change over time. Data <i>volumes</i> of saturation and pressure change values as a function of subsurface location will reduce the uncertainties in 4D interpretations.

00044D seismic is impacted by fluid movement and pressure changes. Production of hydrocarbons resulting in fluid saturation changes in the subsurface will change the acoustic velocity in those regions. Traveling at a different speed (compared to an earlier survey), the reflected seismic wave will arrive at surface detectors sooner or later than in the base survey. This will not only change the strength (amplitude) of the seismic signal, but also shift the apparent depth of reflecting interfaces in resulting seismic images. Moreover, unless pressure is maintained in the reservoir by fluid injection, the reflectors themselves can move downward, which is called subsidence or compaction. However, to quantify the fluid movement and pressure change from the amplitude change of seismic signal, one has to take the rock frame into consideration. Unfortunately, uncertainties in rock parameter estimations (i.e. shale volume vsh and porosity Φ) from seismic can be overwhelming relative to fluid and pressure change signals.

0005Most current technologies for deriving saturation and pressure change use transformations of variations of 4D AVO attributes (quantities calculated from seismic amplitude vs. source-receiver separation, or <i>offset,</i> data), notably A+B and A-B or I<i>_p</i> and I<i>_S</i>. (Data for amplitude vs. offset may be subjected to a best linear fit, i.e. <i>y</i> = A+Bx, in which case the parameters A and B are AVO attributes; I<i>_P</i> and Is represent acoustic impedance, i.e. the products of density and acoustic wave propagation velocities for the P-wave and <i>S</i>-wave.) Using different forms of approximations of the reflectivity equation (see <nplcit id="ncit0001" npl-type="s"><text>Aki and Richards, Quantitative Seismology 123-188 (1980</text></nplcit>)) and a rock physics model, a set of coefficients can be estimated to make a combination of near and far difference amplitudes to infer the saturation and pressure change. See, for example, <nplcit id="ncit0002" npl-type="s"><text>Landro, "Discrimination between pressure and fluid saturation changes from time-lapse seismic data," Geophysics 66, 836-844 (2001</text></nplcit>). This type of formulation works well when physical properties (i.e. vsh and phi) within the reservoir are relatively constant. For a multi-cycle reservoir, side-lobe energy (caused by an input signal such as that from vibroseis for land or airgun for marine being bandwidth limited) may generate apparent difference events that can appear as real reservoir differences. This side-lobe energy complicates interpretation of multi-cycle reservoirs where there is interference between reflectors.

0006Map-based calibration of production data to 4D seismic data has been seen at SEG conventions (<nplcit id="ncit0003" npl-type="s"><text>Floricich et al., "An engineering-driven approach for separating pressure and saturation using 4D seismic: application to a Jurassic reservoir in the UK North Sea," Expanded Abstracts: 75th Annual Meeting of the SEG (2005</text></nplcit>)). Floricich uses production data (pressure and saturation measurements) at well locations and calibrates each quantity with 4D seismic attributes to derive a map. This method is purely statistical. It does not deal with the vertical distribution of saturation and pressure change.

0007Time shifts are commonly used for detecting reservoir compaction. The term <i>time shift</i>, or <i>time lag</i>, refers to comparison of time-lapse seismic data and determining the arrival time correction needed to align the subsurface structure in the later seismic data set with that of the earlier data. See, for example, <nplcit id="ncit0004" npl-type="s"><text>Hudson et al., "Genesis field, Gulf of Mexico, 4-D project status and preliminary lookback," 75th Annual Meeting of the SEG (2005</text></nplcit>); and <nplcit id="ncit0005" npl-type="s"><text>Hatchell et al., "Measuring reservoir compaction using time-lapse timeshifts, Expanded Abstracts," 75th Annual Meeting of the SEG (2005</text></nplcit>).

0008Production will cause the pressure to decrease within the reservoir. If the pressure of the reservoir is not well maintained, compaction of the reservoir will occur, especially for younger rocks. This compaction will most likely be coupled with subsidence of the overburden and overburden dilation. A time-shift data volume can be used to quantify such effects by looking at the time shifts at different time/depth levels.

0009Most recently, at the 2006 SEG meeting, Rickett et al. (2006) and Janssen et al. (2006) used time shifts to estimate the strain caused by production. (<nplcit id="ncit0006" npl-type="s"><text>Rickett et al., "Compacting and 4D time strain at the Genesis Field</text></nplcit>," and <nplcit id="ncit0007" npl-type="s"><text>Janssen et al., "Measuring velocity sensitivity to production-induced strain at the Ekofisk Field using time-lapse time-shifts and compaction logs," both papers in Expanded Abstracts: 76th Annual Meeting of the SEG (2006</text></nplcit>)) Jansen showed that taking the first derivative of the time shifts enables interpreters to interpret time shifts in a manner similar to interpreting 4D difference volumes. Veire et al. disclose a stochastic model for estimation of pressure and saturation changes from time-lapse seismic AVO data within a Bayesian framework. ("<nplcit id="ncit0008" npl-type="s"><text>Stochastic Inversion of Pressure and Saturation Changes From Time-Lapse AVO Data," Geophysics 71, C81 - C92 (2006</text></nplcit>)) XP_002445689 discloses seismic imaging of reservoir flow properties. <patcit id="pcit0002" dnum="US5798982A"><text>US 5,798,982</text></patcit> discloses a method for inverting reflection trace data from 3D and 4D seismic surveys. XP_00124364 discloses direct reservoir parameter estimation using joint inversion of marine seismic AVA and CSEM data.

SUMMARY OF THE INVENTION

0010According to the invention, there is provided a method as defined in any one of the accompanying claims. In one embodiment, the invention is a method for inversion of time-lapse (4D) seismic data (a base survey and a later monitor survey of the same subsurface formation) to estimate the change in fluid saturation and reservoir pressure from the time of the base survey to the time of the monitor survey, comprising: <ol id="ol0001" compact="compact" ol-style=""><li>(a) inverting the 4D seismic data at at least one well location to directly obtain formation pressure change and fluid saturation change for at least one formation fluid, constraining the inversion with well log data; and</li><li>(b) using the saturation and pressure change from at least one well location to obtain saturation and pressure change away from well locations, thereby generating a well-tied fluid saturation and pressure change model predicting change of fluid saturation from the time of the base survey to the time of the monitor survey.</li></ol>

BRIEF DESCRIPTION OF THE DRAWINGS

0011The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which: <ul id="ul0001" list-style="none"><li><figref idref="f0001">Fig. 1</figref> is a workflow chart for one embodiment of the present inventive method;</li><li><figref idref="f0002">Fig. 2</figref> shows inverted (using the present inventive method) saturation and pressure curves for monitor survey data (solid lines) versus corresponding curves for base survey data (dashed lines), wherein the first track is water saturation, the second track is oil saturation, the third track is gas saturation, and the fourth track is reservoir pressure;</li><li><figref idref="f0003">Fig. 3</figref> shows synthetic 4D seismic difference (solid lines), using a model based on inversion results of the present inventive method, versus actual 4D differences (dashed lines), wherein the first track is near angle stack and the second track is far angle stack; and</li><li><figref idref="f0004">Fig. 4</figref> shows (black and white reproduction of a color display) inverted pressure change at the same reservoir represented in <figref idref="f0002 f0003">Figs. 2-3</figref>, the inversion by the present inventive method.</li></ul>

0012The invention will be described in connection with its preferred 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 claims.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

0013The present invention is a method for inferring (i.e., solving an inverse problem) the saturation and pressure change of a reservoir by combining the information from 4D seismic and time lag data volumes derived from the 4D seismic, well logs, and reservoir simulation results (when simulator results are available) and featuring one or more 4D well ties for a quantitative 4D interpretation. The inventive method is model-based incorporating rock physics at well locations, and is statistical-based away from wells. The method thus allows integration of rock physics model and reservoir simulation and honors 4D seismic change. <figref idref="f0001">Figure 1</figref> is a flow chart showing major steps in the workflow in one embodiment of the present invention. In step <b>1</b>, well log data are used to correlate (or <i>tie)</i> events indicated in data from a nearby seismic line with subsurface information provided by the well log. This step is particularly useful, for example, where subsurface layers dip rather than lie parallel to the surface. At step <b>2</b>, rock physics analysis provides values derived from well logs of petrophysical parameters needed for the inversion process - primarily the elastic parameters (<i>P</i>-wave velocity, shear-wave velocity, and density) in relationship to the shale volume fraction (<i>vshale</i> or <i>vsh)</i> and the porosity, with the values valid at least in the vicinity of the well. In step <b>3</b>, the monitor survey data are aligned with base seismic data by apply a volume of time correction factors (i.e. time lag or time shift). Steps <b>1, 2</b> and <b>3</b> are standard steps, well known in the art.

0014At step <b>4</b>, the number of saturation variables is reduced to simplify the inversion step to come. Pressure change (Δ<i>P</i>) is simply the result of subtraction of monitor reservoir pressure from base reservoir pressure. However, saturation change is a bit more complicated. Saturation includes three separate saturation measurements: gas saturation (<i>S_g</i>), oil saturation (<i>S_o</i>) and water saturation (<i>S_w</i>). These three quantities always satisfy the condition of <i>S_g</i>+<i>S_o</i>+<i>S_w</i> = 1. In other words, these quantities are interlinked. To infer saturation change from 4D seismic in a preferred embodiment of the present invention, a single quantity is constructed to represent the changes in all three saturations. It is known that oil saturation rarely increases and that an increase of gas saturation causes a decrease in <i>P</i>-impedance, whereas an increase in water saturation causes an increase in <i>P</i>-impedance. (A <i>P</i>-wave is the primary seismic wave in seismic exploration, and is defined by particle motion in the direction of wave propagation.) Accordingly, the saturation change <i>δ-Sat</i> may be approximated by: <maths id="math0001" num=""><math display="block"><mrow><mi>δ</mi><mo>-</mo><mi mathvariant="italic">Sat</mi><mo>=</mo><mo>-</mo><mi>δ</mi><mo>⁢</mo><msub><mi>S</mi><mi>g</mi></msub><mo>+</mo><mi>d</mi><mo>*</mo><mi>δ</mi><mo>⁢</mo><msub><mi>S</mi><mi>w</mi></msub></mrow></math><img file="EP2153246B1_D0001.tif" /></maths> where <i>d</i> is a correction factor determined by the relative impact of gas versus water on <i>P-</i>impedance, which can be determined through the exercise of Gasmann fluid substitution of oil with brine and gas, separately, and * denotes the multiplication operation. When <i>δ - Sat</i> is negative, this corresponds to an oil saturation decrease and a gas saturation increase; and when <i>δ</i> - <i>Sat</i> is positive, this corresponds to an oil or gas saturation decrease and a water saturation increase.

0015Current mature seismic calibration technology, either volume-based or map-based, deals with one property at a time. The <i>δ - sat</i> of the present invention combines the internal dependency of three saturation changes (gas, oil and water), which makes the seismic calibration for saturation change practical.

0016At step <b>5</b>, model-based inversion is performed at well locations. In other words, saturation and pressure curves are derived by model based-inversion to improve the 4D well tie. The inversion is to find best compromise values for a Δ<i>sat</i> and a Δ<i>P</i> data volume for fitting to 4D difference data and constrained by simulation results (if available) for Δ<i>sat</i> and Δ<i>P</i> based on well log information and a rock physics model. This step is explained in more detail next. One difference between the present inventive method and a published method such as <patcit id="pcit0003" dnum="US5798982A"><text>U.S. Patent No. 5,798,982 to He et al.</text></patcit> is that while He teaches inversion of 4-D seismic data, He inverts to obtain acoustic impedance and then derives oil saturation from the impedance data. He does not directly obtain saturation change from inversion of 4-D seismic data as does the present inventive method.

0017There are several advantages in conducting fluid- and pressure-change inversion at well locations: (a) rock physical properties, as well as elastic properties are best known at well locations; (b) production data and pressure measurements are available only at well locations; and (c) saturation and pressure changes from reservoir simulation are most reliable at well locations. The aforementioned paper by Veire et al. does not teach to conduct inversion at well locations. Further, Veire uses rock physics relationships to develop a stochastic model rather than the present invention's approach of establishing the statistical relationship by directly correlating the inverted saturation and pressure change results with the 4-D seismic data at well locations.

0018In one embodiment of the invention, the inversion scheme involves optimizing the following objective function: <maths id="math0002" num="(1)"><math display="block"><mrow><mi>f</mi><mo>=</mo><mstyle displaystyle="false"><mrow><mstyle displaystyle="true"><mrow><munder><mo>∑</mo><mrow><mi>i</mi></mrow></munder></mrow></mstyle><msup><mfenced separators=""><mi>D</mi><mfenced><msub><mi>m</mi><mi>i</mi></msub></mfenced><mo>-</mo><msub><mi>d</mi><mi>i</mi></msub></mfenced><mrow><mn>2</mn></mrow></msup></mrow></mstyle><mo>+</mo><msub><mrow><mi>λ</mi></mrow><mrow><mn>1</mn></mrow></msub><mo>*</mo><mstyle displaystyle="false"><mrow><mstyle displaystyle="true"><mrow><munder><mo>∑</mo><mrow><mi>i</mi></mrow></munder></mrow></mstyle><msup><mfenced separators=""><msub><mi>m</mi><mi>i</mi></msub><mo>-</mo><msub><mrow><mi>m</mi></mrow><mrow><mi>i</mi><mo>-</mo><mn>1</mn></mrow></msub></mfenced><mrow><mn>2</mn></mrow></msup></mrow></mstyle><mo>+</mo><msub><mrow><mi>λ</mi></mrow><mrow><mn>2</mn></mrow></msub><mo>*</mo><mstyle displaystyle="false"><mrow><mstyle displaystyle="true"><mrow><munder><mo>∑</mo><mrow><mi>i</mi></mrow></munder></mrow></mstyle><msup><mfenced separators=""><msub><mi>m</mi><mi>i</mi></msub><mo>-</mo><msub><mi>r</mi><mi>i</mi></msub></mfenced><mrow><mn>2</mn></mrow></msup></mrow></mstyle></mrow></math><img file="EP2153246B1_D0002.tif" /></maths> where <i>m_i</i> = model saturation and pressure change, <i>d</i> = 4D data, and <i>r_i</i> = reservoir simulation saturation and pressure changes, and the index i is over each discrete depth interval into which the well has been divided for numerical computation purposes. <i>D(m)</i> is the 4D seismic data calculated from the model saturation change. (The term <i>4D data</i> means a 3D data volume of differences in seismic amplitude or some seismic attribute between surveys conducted at different times, time being the 4^th dimension in this terminology.) The λ<i>_i</i> are scaling parameters to balance the terms in the objective function.

0019The first term is a fit to 4D data; the second term controls the smoothness of inverted results; and the third term is a fit to reservoir simulation results. When there is no reservoir simulation available, the following objective function may be used: <maths id="math0003" num="(2)"><math display="block"><mrow><mi>f</mi><mo>=</mo><mstyle displaystyle="false"><mrow><mstyle displaystyle="true"><mrow><munder><mrow><mo>∑</mo></mrow><mrow><mi>i</mi></mrow></munder></mrow></mstyle><msup><mfenced separators=""><mi>D</mi><mfenced><msub><mi>m</mi><mi>i</mi></msub></mfenced><mo>-</mo><msub><mi>d</mi><mi>i</mi></msub></mfenced><mrow><mn>2</mn></mrow></msup></mrow></mstyle><mo>+</mo><msub><mrow><mi>λ</mi></mrow><mrow><mn>1</mn></mrow></msub><mo>*</mo><mstyle displaystyle="false"><mrow><mstyle displaystyle="true"><mrow><munder><mrow><mo>∑</mo></mrow><mrow><mi>i</mi></mrow></munder></mrow></mstyle><msup><mfenced separators=""><msub><mi>m</mi><mi>i</mi></msub><mo>-</mo><msub><mrow><mi>m</mi></mrow><mrow><mi>i</mi><mo>-</mo><mn>1</mn></mrow></msub></mfenced><mrow><mn>2</mn></mrow></msup></mrow></mstyle><mo>+</mo><msub><mrow><mi>λ</mi></mrow><mrow><mn>2</mn></mrow></msub><mo>*</mo><msup><mrow><mfenced open="[" close="]" separators=""><mstyle displaystyle="false"><mrow><mstyle displaystyle="true"><mrow><munder><mo>∑</mo><mrow><mi>i</mi></mrow></munder></mrow></mstyle><mrow><mfenced separators=""><mi>m</mi><mo>*</mo><mfenced separators=""><mn>1</mn><mo>-</mo><mi mathvariant="italic">vsh</mi></mfenced><mo>*</mo><mi>φ</mi></mfenced></mrow></mrow></mstyle><mo>-</mo><mi mathvariant="italic">δv</mi></mfenced></mrow><mrow><mn>2</mn></mrow></msup></mrow></math><img file="EP2153246B1_D0003.tif" /></maths> where δ<i>v</i> is the ratio of produced gas-to-oil or water-to-oil. The third term is to match the produced portions of liquids (water, gas and water). Skilled readers will be able to devise other choices for the objective functions above.

0020Next, a way to represent the objective function as a set of linear equations will be given. If it is assumed that <i>P</i>-wave velocity (<i>V_p</i>), <i>S</i>-wave velocity (<i>V_S</i>), and density (<i>ρ</i>) are known values of those elastic parameters before production began, the reflectivity as a function of angle (θ) for the base seismic survey can be written as follows (see <nplcit id="ncit0009" npl-type="s"><text>Stolt and Weglein, "Migration and inversion of seismic data," Geophysics 50, 2458-2472 (1985</text></nplcit>) and the previously cited Aki &amp; Richards 1980 reference): <maths id="math0004" num="(3)"><math display="block"><mrow><mtable columnalign="left"><mtr><mtd><mi mathvariant="italic">Ref</mi><mfenced separators=","><mi>t</mi><mi>θ</mi></mfenced></mtd><mtd><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mn>2</mn></mrow></mfrac><mo>⁢</mo><mfenced separators=""><mn>1</mn><mo>+</mo><msup><mrow><mi>tan</mi></mrow><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><mi>θ</mi></mfenced><mrow><mfrac><mrow><mo>∂</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac></mrow><mi>ln</mi><mrow><mi mathvariant="italic">Vp</mi></mrow><mfenced><mi>t</mi></mfenced><mo>-</mo><mn>4</mn><mo>⁢</mo><mfrac><mrow><msup><mrow><mover><mrow><msub><mi>V</mi><mi>s</mi></msub></mrow><mrow><mo>‾</mo></mrow></mover></mrow><mrow><mn>2</mn></mrow></msup></mrow><mrow><msup><mrow><mover><mrow><msub><mi>V</mi><mi>p</mi></msub></mrow><mrow><mo>‾</mo></mrow></mover></mrow><mrow><mn>2</mn></mrow></msup></mrow></mfrac><mo>⁢</mo><msup><mrow><mi>sin</mi></mrow><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><mi>θ</mi><mrow><mfrac><mrow><mo>∂</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac></mrow><mi>ln</mi><mrow><mi mathvariant="italic">Vs</mi></mrow><mfenced><mi>t</mi></mfenced></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mo>+</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mn>2</mn></mrow></mfrac><mo>⁢</mo><mfenced separators=""><mn>1</mn><mo>-</mo><mn>4</mn><mo>⁢</mo><mfrac><mrow><msup><mrow><mover><mrow><msub><mi>V</mi><mi>s</mi></msub></mrow><mrow><mo>‾</mo></mrow></mover></mrow><mrow><mn>2</mn></mrow></msup></mrow><mrow><msup><mrow><mover><mrow><msub><mi>V</mi><mi>p</mi></msub></mrow><mrow><mo>‾</mo></mrow></mover></mrow><mrow><mn>2</mn></mrow></msup></mrow></mfrac><mo>⁢</mo><msup><mrow><mi>sin</mi></mrow><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><mi>θ</mi></mfenced><mo>⁢</mo><mfrac><mrow><mo>∂</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>⁢</mo><mi>ln</mi><mrow><mi>ρ</mi></mrow><mfenced><mi>t</mi></mfenced></mtd></mtr></mtable></mrow></math><img file="EP2153246B1_D0004.tif" /></maths>

0021Seismic data are obtained by convolving the reflectivity with a wavelet. This can be written as a set of matrix equations: <maths id="math0005" num="(4)"><math display="block"><mrow><mfenced open="[" close="]"><mtable><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mi mathvariant="italic">Wavelet</mi></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr></mtable></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mi mathvariant="italic">Aki</mi><mo>&amp;</mo></mtd></mtr><mtr><mtd><mi mathvariant="italic">Richards</mi></mtd></mtr></mtable></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mi mathvariant="italic">Derivative</mi></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr></mtable></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">Vp</mi></mrow></mtd></mtr><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">Vs</mi></mrow></mtd></mtr><mtr><mtd><mi>ln</mi><mspace width="1em" /><mrow><mi mathvariant="italic">ρ</mi></mrow></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open="[" close="]"><mi mathvariant="italic">base</mi></mfenced></mrow></math><img file="EP2153246B1_D0005.tif" /></maths> that can be rewritten simply as a single matrix "G": <maths id="math0006" num="(5)"><math display="block"><mrow><mfenced open="[" close="]" separators=""><mspace width="1em" /><mtable><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mi>G</mi></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr></mtable><mspace width="1em" /></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">Vp</mi></mrow></mtd></mtr><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">Vs</mi></mrow></mtd></mtr><mtr><mtd><mi>ln</mi><mrow><mspace width="1em" /></mrow><mrow><mi mathvariant="italic">ρ</mi></mrow></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open="[" close="]"><mi mathvariant="italic">base</mi></mfenced><mn>.</mn></mrow></math><img file="EP2153246B1_D0006.tif" /></maths>

0022After production, for an incremental change of pressure and saturation, the effect of saturation change on shear wave and effect of pressure on density may both be considered negligible, resulting in: <maths id="math0007" num="(6)"><math display="block"><mrow><mtable columnalign="left"><mtr><mtd><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi mathvariant="italic">Vp</mi></mrow><mo>≈</mo><mi mathvariant="italic">Ap</mi><mo>*</mo><mi mathvariant="italic">dsat</mi><mo>+</mo><mi mathvariant="italic">Bp</mi><mo>*</mo><mi mathvariant="italic">dp</mi></mtd></mtr><mtr><mtd><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi mathvariant="italic">Vs</mi></mrow><mo>≈</mo><mn>0</mn><mo>*</mo><mi mathvariant="italic">dsat</mi><mo>+</mo><mi mathvariant="italic">Bs</mi><mo>*</mo><mi mathvariant="italic">dp</mi></mtd></mtr><mtr><mtd><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi mathvariant="italic">ρ</mi></mrow><mo>≈</mo><msub><mi mathvariant="italic">A</mi><mi>ρ</mi></msub><mo>*</mo><mi mathvariant="italic">dsat</mi><mo>+</mo><mn>0</mn><mo>*</mo><mi mathvariant="italic">dp</mi></mtd></mtr></mtable></mrow></math><img file="EP2153246B1_D0007.tif" /></maths>

0023Gassmann fluid substitution and an empirical relationship may be used to derive the coefficients <i>A_P, B_P, B_S</i> and <i>A_ρ</i> (see Appendix). For monitor seismic data: <maths id="math0008" num="(7)"><math display="block"><mrow><mrow><mfenced open="[" close="]" separators=""><mspace width="1em" /><mtable><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mi>G</mi></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr></mtable><mspace width="1em" /></mfenced></mrow><mfenced open="[" close="]"><mtable columnalign="left"><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">Vp</mi></mrow><mo>+</mo><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi mathvariant="italic">Vp</mi></mrow></mtd></mtr><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">Vs</mi></mrow><mo>+</mo><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi mathvariant="italic">Vs</mi></mrow></mtd></mtr><mtr><mtd><mi>ln</mi><mrow><mi mathvariant="italic">ρ</mi></mrow><mo>+</mo><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi mathvariant="italic">ρ</mi></mrow></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open="[" close="]"><mi mathvariant="italic">Monitor</mi></mfenced><mn>.</mn></mrow></math><img file="EP2153246B1_D0008.tif" /></maths>

0024Therefore <maths id="math0009" num="(8)"><math display="block"><mrow><mrow><mfenced open="[" close="]"><mtable><mtr><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mi>G</mi></mtd><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd></mtr></mtable></mfenced></mrow><mfenced open="[" close="]"><mtable><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mi>A</mi></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd></mtr></mtable></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable columnalign="left"><mtr><mtd><mi mathvariant="italic">dsat</mi></mtd></mtr><mtr><mtd><mi mathvariant="italic">dp</mi></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open="[" close="]" separators=""><mi mathvariant="italic">Monitor</mi><mo>-</mo><mi mathvariant="italic">base</mi></mfenced><mn>.</mn></mrow></math><img file="EP2153246B1_D0009.tif" /></maths> Where <maths id="math0010" num=""><math display="inline"><mrow><mi>A</mi><mo>=</mo><mfenced open="[" close="]"><mtable columnalign="left"><mtr><mtd><mi>A</mi><mrow><mi>p</mi></mrow></mtd><mtd><mi>B</mi><mrow><mi>p</mi></mrow></mtd></mtr><mtr><mtd><mn>0</mn></mtd><mtd><mi>B</mi><mrow><mi>s</mi></mrow></mtd></mtr><mtr><mtd><mi>A</mi><mrow><mi>ρ</mi></mrow></mtd><mtd><mn>0</mn></mtd></mtr></mtable></mfenced><mn>.</mn></mrow></math><img file="EP2153246B1_D0010.tif" /></maths>

0025As will be well understood by practitioners in the field, the inversion uses an iterative process with each inversion step to solve for incremental improvement of goodness of fit.

0026It was found that it is typically most stable to directly apply the present inventive method to the 4D difference data rather than separately inverting monitor data and base data and then comparing. For objective function based inversion, a damped least-squares formalism may be used to solve the linearized set of equations. This involves appending two sets of damping matrices below the <i>G*A</i> matrix in equation (<b>8</b>). The first set of damping parameters control smoothness of the solution, which are differential operators across the diagonal. The second set of damping parameters are for fitting to the reservoir simulation. The damping matrix and associated portion of the data vector are multiplied by coefficients that are determined by a series of tests in which the coefficients are varied. The optimal damping parameters are chosen and the resulting matrix has the following form: <maths id="math0011" num="(9)"><math display="block"><mrow><mfenced open="[" close="]"><mtable><mtr><mtd><mi>G</mi><mo>*</mo><mi>A</mi></mtd></mtr><mtr><mtd><msub><mrow><mi>λ</mi></mrow><mrow><mn>1</mn></mrow></msub><mo>*</mo><mi>D</mi></mtd></mtr><mtr><mtd><msub><mrow><mi>λ</mi></mrow><mrow><mn>2</mn><mo>⁢</mo><mi>s</mi></mrow></msub><mo>*</mo><msub><mrow><mi>I</mi></mrow><mrow><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>λ</mi></mrow><mrow><mn>2</mn><mo>⁢</mo><mi>p</mi></mrow></msub><mo>*</mo><msub><mrow><mi>I</mi></mrow><mrow><mn>2</mn></mrow></msub></mtd></mtr></mtable></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable columnalign="left"><mtr><mtd><mi mathvariant="italic">dsat</mi></mtd></mtr><mtr><mtd><mi mathvariant="italic">dp</mi></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open="[" close="]"><mtable columnalign="left"><mtr><mtd><mi mathvariant="italic">difference</mi></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mi mathvariant="italic">dsat_res</mi></mtd></mtr><mtr><mtd><mi mathvariant="italic">dp_res</mi></mtd></mtr></mtable></mfenced></mrow></math><img file="EP2153246B1_D0011.tif" /></maths> with <maths id="math0012" num=""><math display="inline"><mrow><mi>D</mi><mo>=</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mo>…</mo></mtd><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mo>-</mo><mn>1</mn></mtd><mtd><mn>1</mn></mtd><mtd><mspace width="1em" /></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd><mtd><mo>-</mo><mn>1</mn></mtd><mtd><mn>1</mn></mtd></mtr><mtr><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd><mtd><mspace width="1em" /></mtd><mtd><mo>…</mo></mtd></mtr></mtable></mfenced><mo>,</mo></mrow></math><img file="EP2153246B1_D0012.tif" /></maths><i>I</i>₁ = [<i>I</i> 0], <i>I</i>₂ = [0 <i>I</i>] and <i>I</i> is an identity matrix.

0027With appropriate choices of λ, this particular parameterization produces a solution that improves the fit to 4D data by honoring the reservoir simulation and a smooth change of saturation and pressure.

0028For the subsequent iteration of the inversion, the following form may be used: <maths id="math0013" num="(10)"><math display="block"><mrow><mfenced open="[" close="]"><mtable><mtr><mtd><mi>G</mi><mo>*</mo><mi>A</mi></mtd></mtr><mtr><mtd><msub><mrow><mi>λ</mi></mrow><mrow><mn>1</mn></mrow></msub><mo>*</mo><mi>D</mi></mtd></mtr><mtr><mtd><msub><mrow><mi>λ</mi></mrow><mrow><mn>2</mn><mo>⁢</mo><mi>s</mi></mrow></msub><mo>*</mo><msub><mrow><mi>I</mi></mrow><mrow><mn>1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>λ</mi></mrow><mrow><mn>2</mn><mo>⁢</mo><mi>p</mi></mrow></msub><mo>*</mo><msub><mrow><mi>I</mi></mrow><mrow><mn>2</mn></mrow></msub></mtd></mtr></mtable></mfenced><mo>⁢</mo><mfenced open="[" close="]"><mtable><mtr><mtd><mi mathvariant="italic">ddsat</mi></mtd></mtr><mtr><mtd><mi mathvariant="italic">ddp</mi></mtd></mtr></mtable></mfenced><mo>=</mo><mfenced open="[" close="]"><mtable columnalign="left"><mtr><mtd><mi mathvariant="italic">seismic diff</mi><mo>-</mo><mi mathvariant="italic">modeled diff</mi></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable></mfenced></mrow></math><img file="EP2153246B1_D0013.tif" /></maths> where <i>ddsat</i> and <i>ddp</i> are incremental updates to <i>dsat</i> and <i>dp.</i> Three to four iterations are in general enough to obtain a converged solution.

0029After saturation and fluid change inversions are completed at well locations, pressure and saturation change away from well locations can be obtained (step <b>6</b>) through either one of at least two approaches. The first approach is using the rock physics and inversion parameters derived at well locations to conduct fluid and pressure change inversions at pseudo-well locations. This approach requires a high confidence about the rock properties (i.e. <i>vsh</i> and porosity) inverted from seismic data at pseudo-well locations. The second approach is to use statistical calibration, such as multivariate linear regression or neural networks, to extrapolate away from well locations. Inverted saturation and pressure changes at well locations are calibrated to 4D seismic differences from different offsets. This approach requires wells used for calibration covering the rock property variations at reservoir.

0030In step 7, a low-frequency saturation change model is generated from time lag data volumes. (Before the time lag volumes are generated, a standard 4D seismic step is performed to align the time scales for the base seismic data and the monitor seismic data. This is shown in <figref idref="f0001">Fig. 1</figref> as step 3. The time alignment is also required for step 5.)

0031In one embodiment of step 7, let <i>T_lag</i> be the time lag required to shift the monitor survey to align it with the base survey at time <i>t</i>. It was found that the derivative of the time lag data volume is an approximation of the instant relative velocity change (<i>δv</i>/<i>v</i>) between monitor and base when there is no significant compaction of at reservoir: <maths id="math0014" num="(11)"><math display="block"><mrow><mrow><mfrac><mrow><mo>ⅆ</mo><msub><mi>T</mi><mi mathvariant="italic">lag</mi></msub></mrow><mrow><mo>ⅆ</mo><mi>t</mi></mrow></mfrac></mrow><mo>≈</mo><mfrac><mrow><mi mathvariant="italic">δv</mi></mrow><mrow><mn>2</mn><mo>⁢</mo><mi>v</mi></mrow></mfrac></mrow></math><img file="EP2153246B1_D0014.tif" /></maths> where <i>δv = v_m</i> - <i>v_b</i>, <i>v_m</i> = velocity at the time of the monitor survey and <i>v_b</i> = base survey velocity. For <i>P</i>-impedance <i>(I_P): Ip = v*ρ</i>; therefore, <maths id="math0015" num="(12)"><math display="block"><mrow><mfrac><mrow><mi>δ</mi><mo>⁢</mo><msub><mi>I</mi><mi>p</mi></msub></mrow><mrow><msub><mi>I</mi><mi>p</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><mi mathvariant="italic">δv</mi></mrow><mrow><mi>v</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mi mathvariant="italic">δρ</mi></mrow><mrow><mi>ρ</mi></mrow></mfrac></mrow></math><img file="EP2153246B1_D0015.tif" /></maths> On the assumption that the density contribution is small compared to the velocity contribution: <maths id="math0016" num=""><math display="block"><mrow><mi mathvariant="italic">δIp</mi><mo>=</mo><mi mathvariant="italic">Ip</mi><mo>*</mo><mfrac><mrow><mi mathvariant="italic">δv</mi></mrow><mrow><mi>v</mi></mrow></mfrac></mrow></math><img file="EP2153246B1_D0016.tif" /></maths>

0032Substitution of equation 11 into the last equation yields the time-lag based impedance change. <i>δI_P,</i> the impedance change, is in the same domain as the amplitude difference, i.e. both are a measure of reflection strength. When an inversion is not available, integrated absolute seismic trace amplitude <i>(I_AA)</i> may be used as a replacement of <i>I_P.</i> The potential applications of time-lag based impedance are many, including: (a) to provide a low-frequency model for difference-based inversion; (b) to modify the base survey low-frequency model for improved monitor impedance inversion; and (c) to display as a cross-plot with the amplitude-based impedance change for improved 4D reservoir characterizations.

0033The low frequency model may be created by calibrating the low-filtered saturation change logs, which were derived in step <b>5</b> with a similarly low-pass filtered <i>δI_P</i> volume.

0034In step <b>8</b>, the low frequency model may be merged with a band-limited saturation change volume derived from the calibration step to get a high-cut saturation volume; i.e. the two data volumes are added together after making sure there is little overlap of the information in frequency domain. The merged volume improves correspondence to the water/gas movement and accuracy of the saturation change predication, as indicated by the results in Table 1.

0035Inverted saturation and pressure change can be used to more rigorously constrain the reservoir simulation model or as target logs for statistical calibration to extrapolate away from wells, as is done in the example that follows.

Example

0036The present inventive method was applied to a deepwater field where gas and water injection is used for both pressure maintenance and gas storage. Monitoring the gas and water movement is the main objective of 4D surveillance at this field. <figref idref="f0002">Figure 2</figref> shows the inverted saturation and pressure change curves at a well location. The solid lines correspond to the monitor survey, and the dashed lines the base survey. <figref idref="f0003">Figure 3</figref> plots synthetic 4D seismic difference (solid lines) versus actual 4D differences (dashed lines). The first track is a near angle stack and second track is a far angle stack. <figref idref="f0003">Figure 3</figref> shows that the fluid and pressure change inversion of the present invention (solid lines) is able to match actual differences very well for the near angle stack, and less well for the far angle stack (to be expected because of factors such as uncertainties in NMO stretch and more Q attenuation).

0037Using inverted logs at this field as target logs, the commercial software product Emerge (http://www.veritasdgc.com/vhr/bins/index.asp) was used to calibrate 4D seismic to saturation and pressure change. Both near and far difference volumes are among the input volumes. As <figref idref="f0004">Fig. 4</figref> shows for Δ<i>P</i> and Table 1 shows for saturation, the resultant pressure and saturation show good conformance to geological expectations and consistency with reservoir simulations. <figref idref="f0004">Figure 4</figref> is a black and white reproduction of a colored data display where the color scale represents numerical values of Δ<i>P</i>. <tables id="tabl0001" num="0001"><table frame="all"><title>Table 1.</title><tgroup cols="3"><colspec colnum="1" colname="col1" colwidth="68mm" /><colspec colnum="2" colname="col2" colwidth="27mm" /><colspec colnum="3" colname="col3" colwidth="14mm" /><thead><row><entry valign="top">Volumes</entry><entry valign="top">Correlation (R)</entry><entry valign="top">Error</entry></row></thead><tbody><row><entry>Far angle stack</entry><entry>0.20</entry><entry>0.102</entry></row><row><entry>Bandlimited Saturation volume from inversion</entry><entry>0.32</entry><entry>0.096</entry></row><row><entry>Merged Saturation Volume</entry><entry>0.78</entry><entry>0.066</entry></row></tbody></tgroup></table></tables>

0038The fluid saturation and pressure change versus depth results can be used to improve the history match in reservoir simulation, i.e. to provide added constraint on the reservoir simulation model.

0039The foregoing 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. 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 recognize that at least some of the steps of the invention are performed on or with the aid of a computer, i.e. the invention is computer implemented in its preferred embodiments. In such instances, the resulting pressure and saturation changes may be either downloaded or saved to computer memory.

Appendix

0040To derive the coefficients in equation (6), it is assumed that pressure change and saturation change affect the elastic parameters (<i>V_P, V_S</i> and <i>ρ</i>) independently. The coefficients for saturation and pressure are derived using separate equations:

<u>Saturation change</u>

0041The effect of saturation change on shear velocity is ignored. For P-velocity: <maths id="math0017" num="(A1)"><math display="block"><mrow><msub><mi>V</mi><mi>p</mi></msub><mo>=</mo><msqrt><mrow><mfrac><mrow><mi>K</mi><mo>+</mo><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><mn>3</mn></mrow></msub></mrow><mprescripts /><none /><mrow><mn>4</mn></mrow></mmultiscripts><mi>μ</mi></mrow><mrow><mi>ρ</mi></mrow></mfrac></mrow></msqrt></mrow></math><img file="EP2153246B1_D0017.tif" /></maths> Only bulk modulus (<i>K</i>) and density (<i>ρ</i>) vary when saturation changes, therefore <maths id="math0018" num="(A2)"><math display="block"><mrow><mtable><mtr><mtd><mi>d</mi><mrow><mi>ln</mi></mrow><msub><mi>V</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><msub><mi>V</mi><mi>p</mi></msub></mrow></mfrac><mo>⁢</mo><mfenced separators=""><mrow><mfrac><mrow><mo>ⅆ</mo><msub><mi>V</mi><mi>p</mi></msub></mrow><mrow><mo>ⅆ</mo><mi>K</mi></mrow></mfrac></mrow><mi mathvariant="italic">dK</mi><mo>+</mo><mrow><mfrac><mrow><mo>ⅆ</mo><msub><mi>V</mi><mi>p</mi></msub></mrow><mrow><mo>ⅆ</mo><mi>ρ</mi></mrow></mfrac></mrow><mi mathvariant="italic">dρ</mi></mfenced></mtd></mtr><mtr><mtd><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mn>2</mn><mo>⁢</mo><msup><mrow><msub><mi>V</mi><mi>p</mi></msub></mrow><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><mi>ρ</mi></mrow></mfrac><mo>⁢</mo><mfenced separators=""><mi mathvariant="italic">dK</mi><mo>-</mo><mi>d</mi><mrow><mi>ln</mi></mrow><mi>ρ</mi></mfenced></mtd></mtr></mtable></mrow></math><img file="EP2153246B1_D0018.tif" /></maths> and <maths id="math0019" num="(A3)"><math display="block"><mrow><mi>d</mi><mrow><mi>ln</mi></mrow><mfenced><mi>ρ</mi></mfenced><mo>=</mo><mn>1</mn><mo>/</mo><mi mathvariant="italic">ρdρ</mi></mrow></math><img file="EP2153246B1_D0019.tif" /></maths>

0042From equation (A3) and the following two equations: <maths id="math0020" num=""><math display="block"><mrow><mi>ρ</mi><mo>=</mo><mfenced separators=""><mn>1</mn><mo>-</mo><mi>φ</mi></mfenced><mo>⁢</mo><msub><mi>ρ</mi><mi>s</mi></msub><mo>+</mo><mi>φ</mi><mo>⁢</mo><msub><mi>ρ</mi><mi>f</mi></msub><mspace width="1em" /><mrow><mi mathvariant="italic">and</mi></mrow><mspace width="1em" /><msub><mi>ρ</mi><mi>f</mi></msub><mo>=</mo><msub><mi>ρ</mi><mi>o</mi></msub><mo>⁢</mo><msub><mi>S</mi><mi>o</mi></msub><mo>+</mo><msub><mi>ρ</mi><mi>g</mi></msub><mo>⁢</mo><msub><mi>S</mi><mi>g</mi></msub></mrow></math><img file="EP2153246B1_D0020.tif" /></maths> where <i>ρ_s</i> is the density for solid rock, <i>ρ_f</i> is the density for mixed fluid, <i>ρ_o</i> is the density of in-situ oil and <i>ρ_g</i> is the density of in-situ gas, one then has: <maths id="math0021" num="(A4)"><math display="block"><mrow><mi>d</mi><mrow><mi>ln</mi></mrow><mrow><mi>ρ</mi></mrow><mo>=</mo><mfrac><mrow><mi>φ</mi><mo>⁢</mo><mfenced separators=""><msub><mi>ρ</mi><mi>o</mi></msub><mo>-</mo><msub><mi>ρ</mi><mi>g</mi></msub></mfenced></mrow><mrow><mi>ρ</mi></mrow></mfrac><mo>⁢</mo><mi mathvariant="italic">dsat</mi></mrow></math><img file="EP2153246B1_D0021.tif" /></maths> Here, only a two-fluid system of gas-oil will be considered because the correction factor d in equation (1) is designed to account for water replacing oil. Hence, <maths id="math0022" num=""><math display="inline"><mrow><msub><mi>A</mi><mi>ρ</mi></msub><mo>=</mo><mfrac><mrow><mi>φ</mi><mo>⁢</mo><mfenced separators=""><msub><mi>ρ</mi><mi>o</mi></msub><mo>-</mo><msub><mi>ρ</mi><mi>g</mi></msub></mfenced></mrow><mrow><mi>ρ</mi></mrow></mfrac></mrow></math><img file="EP2153246B1_D0022.tif" /></maths>

0043To express <i>dK</i> in terms of <i>dsat,</i> one may start with the Gassmann equation: <maths id="math0023" num="(A5)"><math display="block"><mrow><mfrac><mrow><msub><mi>K</mi><mi mathvariant="italic">satA</mi></msub></mrow><mrow><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><msub><mi>K</mi><mi mathvariant="italic">satA</mi></msub></mrow></mfrac><mo>-</mo><mfrac><mrow><msub><mi>K</mi><mi mathvariant="italic">fluidA</mi></msub></mrow><mrow><mi>φ</mi><mo>⁢</mo><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><msub><mi>K</mi><mi mathvariant="italic">fluidA</mi></msub></mfenced></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>K</mi><mi mathvariant="italic">satB</mi></msub></mrow><mrow><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><msub><mi>K</mi><mi mathvariant="italic">satB</mi></msub></mrow></mfrac><mo>-</mo><mfrac><mrow><msub><mi>K</mi><mi mathvariant="italic">fluidB</mi></msub></mrow><mrow><mi>φ</mi><mo>⁢</mo><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><msub><mi>K</mi><mi mathvariant="italic">fluidB</mi></msub></mfenced></mrow></mfrac></mrow></math><img file="EP2153246B1_D0023.tif" /></maths><i>Assume k = K_satA</i> and <i>k_f</i> = <i>K_fluidA</i>; for the rock frame with fluid <i>B,</i> one has <i>K_satB = k+dk</i> and <i>K_fluidB = k</i>_<i>f</i>+<i>dk_f.</i> After plugging these definitions into equation (A5) and neglecting small terms, one gets: <maths id="math0024" num="(A6)"><math display="block"><mrow><mi mathvariant="italic">dk</mi><mo>=</mo><mfrac><mrow><msup><mrow><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><mi>K</mi></mfenced></mrow><mrow><mn>2</mn></mrow></msup></mrow><mrow><mi>φ</mi><mo>⁢</mo><msup><mrow><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><msub><mi>K</mi><mi>f</mi></msub></mfenced></mrow><mrow><mn>2</mn></mrow></msup></mrow></mfrac><mo>⁢</mo><msub><mi mathvariant="italic">dk</mi><mi>f</mi></msub></mrow></math><img file="EP2153246B1_D0024.tif" /></maths>

0044Similar to density, a two-fluid system is used to express <i>dk_f</i> in terms of <i>dsat.</i> For a Reuss average mixing scheme: <maths id="math0025" num="(A7)"><math display="block"><mrow><mfrac><mrow><mn>1</mn></mrow><mrow><msub><mi>K</mi><mi>f</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mrow><msub><mi>S</mi><mi>g</mi></msub></mrow><mrow><msub><mi>K</mi><mi>g</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>S</mi><mi>o</mi></msub></mrow><mrow><msub><mi>K</mi><mi>o</mi></msub></mrow></mfrac></mrow></math><img file="EP2153246B1_D0025.tif" /></maths> One finds: <maths id="math0026" num="(A8)"><math display="block"><mrow><mi>d</mi><mo>⁢</mo><msub><mi>K</mi><mi>f</mi></msub><mo>=</mo><mfrac><mrow><msubsup><mrow><mi>K</mi></mrow><mrow><mi>f</mi></mrow><mrow><mn>2</mn></mrow></msubsup><mo>⁢</mo><mfenced separators=""><msub><mi>K</mi><mi>o</mi></msub><mo>-</mo><msub><mi>K</mi><mi>g</mi></msub></mfenced></mrow><mrow><msub><mi>K</mi><mi>o</mi></msub><mo>⁢</mo><msub><mi>K</mi><mi>g</mi></msub></mrow></mfrac><mo>⁢</mo><mi mathvariant="italic">dsat</mi></mrow></math><img file="EP2153246B1_D0026.tif" /></maths> By substituting equations (A8), (A6) and (A4) into (A2), one gets: <maths id="math0027" num="(A9)"><math display="block"><mrow><mi>d</mi><mrow><mi>ln</mi></mrow><mi mathvariant="italic">Vp</mi><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mn>2</mn><mo>⁢</mo><msubsup><mrow><mi>V</mi></mrow><mrow><mi>p</mi></mrow><mrow><mn>2</mn></mrow></msubsup><mo>⁢</mo><mi>ρ</mi></mrow></mfrac><mo>⁢</mo><mfenced separators=""><mfrac><mrow><msup><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><mi>K</mi></mfenced><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><msubsup><mrow><mi>K</mi></mrow><mrow><mi>f</mi></mrow><mrow><mn>2</mn></mrow></msubsup><mo>⁢</mo><mfenced separators=""><msub><mi>K</mi><mi>o</mi></msub><mo>-</mo><msub><mi>K</mi><mi>g</mi></msub></mfenced></mrow><mrow><mi>φ</mi><mrow><msup><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><msub><mi>K</mi><mi>f</mi></msub></mfenced><mrow><mn>2</mn></mrow></msup></mrow><msub><mi>K</mi><mi>o</mi></msub><mo>⁢</mo><msub><mi>K</mi><mi>g</mi></msub></mrow></mfrac><mo>-</mo><mfrac><mrow><mi>φ</mi><mo>⁢</mo><mfenced separators=""><msub><mi>ρ</mi><mi>o</mi></msub><mo>-</mo><msub><mi>ρ</mi><mi>g</mi></msub></mfenced></mrow><mrow><mi>ρ</mi></mrow></mfrac></mfenced><mo>⁢</mo><mi mathvariant="italic">dsat</mi></mrow></math><img file="EP2153246B1_D0027.tif" /></maths> Hence, <maths id="math0028" num=""><math display="inline"><mrow><msub><mi>A</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><mn>1</mn></mrow><mrow><mn>2</mn><mo>⁢</mo><msubsup><mrow><mi>V</mi></mrow><mrow><mi>p</mi></mrow><mrow><mn>2</mn></mrow></msubsup><mo>⁢</mo><mi>ρ</mi></mrow></mfrac><mo>⁢</mo><mfenced separators=""><mfrac><mrow><msup><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><mi>K</mi></mfenced><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><msubsup><mrow><mi>K</mi></mrow><mrow><mi>f</mi></mrow><mrow><mn>0</mn></mrow></msubsup><mo>⁢</mo><mfenced separators=""><msub><mi>K</mi><mi>o</mi></msub><mo>-</mo><msub><mi>K</mi><mi>g</mi></msub></mfenced></mrow><mrow><mi>φ</mi><mo>⁢</mo><msup><mfenced separators=""><msub><mi>K</mi><mi>s</mi></msub><mo>-</mo><mi>K</mi></mfenced><mrow><mn>2</mn></mrow></msup><mo>⁢</mo><msub><mi>K</mi><mi>o</mi></msub><mo>⁢</mo><msub><mi>K</mi><mi>g</mi></msub></mrow></mfrac><mo>-</mo><mfrac><mrow><mi>φ</mi><mo>⁢</mo><mfenced separators=""><msub><mi>ρ</mi><mi>o</mi></msub><mo>-</mo><msub><mi>ρ</mi><mi>g</mi></msub></mfenced></mrow><mrow><mi>ρ</mi></mrow></mfrac></mfenced></mrow></math><img file="EP2153246B1_D0028.tif" /></maths>

<u>Pressure Change</u>

0045The pressure effect on density will be ignored, and empirical relationships will be used for <i>V_P</i> and <i>V_S.</i> For <i>Vp</i> in terms of pressure <i>P,</i><maths id="math0029" num="(A10)"><math display="block"><mrow><mi mathvariant="italic">Vp</mi><mo>=</mo><msub><mi mathvariant="italic">Vp</mi><mi>∞</mi></msub><mo>⁢</mo><mfenced separators=""><mn>1</mn><mo>-</mo><mi>A</mi><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mfenced></mrow></math><img file="EP2153246B1_D0029.tif" /></maths><i>A</i> and <i>P₀</i> are constants, which are generally determined from core measurements. By taking the logarithm and the derivative relative to <i>P</i>, one gets <maths id="math0030" num="(A11)"><math display="block"><mrow><mi mathvariant="italic">d</mi><mrow><mi>ln</mi></mrow><mi mathvariant="italic">Vp</mi><mo>=</mo><mfrac><mrow><mi>A</mi><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mrow><mrow><msub><mrow><mi>p</mi></mrow><mrow><mn>0</mn></mrow></msub><mo>⁢</mo><mfenced separators=""><mn>1</mn><mo>-</mo><mi>A</mi><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mfenced></mrow></mfrac><mo>⁢</mo><mi mathvariant="italic">dp</mi></mrow></math><img file="EP2153246B1_D0030.tif" /></maths> Hence, <maths id="math0031" num=""><math display="block"><mrow><msub><mi>B</mi><mi>p</mi></msub><mo>=</mo><mfrac><mrow><mi>A</mi><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mrow><mrow><msub><mrow><mi>p</mi></mrow><mrow><mn>0</mn></mrow></msub><mo>⁢</mo><mfenced separators=""><mn>1</mn><mo>-</mo><mi>A</mi><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mfenced></mrow></mfrac></mrow></math><img file="EP2153246B1_D0031.tif" /></maths>

0046Similarly for shear velocity <i>V_S,</i> one has <maths id="math0032" num=""><math display="block"><mrow><mi mathvariant="italic">Bs</mi><mo>=</mo><mfrac><mrow><msub><mi>A</mi><mi>s</mi></msub><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mi>s</mi><mo>⁢</mo><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mrow><mrow><msub><mrow><mi>p</mi></mrow><mrow><mi>s</mi><mo>⁢</mo><mn>0</mn></mrow></msub><mo>⁢</mo><mfenced separators=""><mn>1</mn><mo>-</mo><msub><mi>A</mi><mi>s</mi></msub><mo>⁢</mo><msup><mrow><mi>e</mi></mrow><mrow><mmultiscripts><mrow><msub><mrow><mo>/</mo></mrow><mrow><msub><mrow><mi>P</mi></mrow><mrow><mi>s</mi><mo>⁢</mo><mn>0</mn></mrow></msub></mrow></msub></mrow><mprescripts /><none /><mrow><mo>-</mo><mi>p</mi></mrow></mmultiscripts></mrow></msup></mfenced></mrow></mfrac><mo>⁢</mo><mi mathvariant="italic">dp</mi></mrow></math><img file="EP2153246B1_D0032.tif" /></maths>