Research on Remaining Oil Characterization in Superheavy Oil Reservoir by Microgravity Exploration

Abstract

Some physical processes such as oil and gas development, metal deposit collection, and groundwater resource migration can cause density changes, for which microgravity monitoring is the most intuitive method to monitor the density change process. Based on the basic principle of microgravity measurement and the idea of multiscale separation, a multiscale, second-order, surface-fitting, residual gravity anomaly extraction method is proposed to separate superimposed microgravity fields. In this method, regional fields of different scales are fitted and calculated successively with the measurement points as the center, so as to separate the gravity anomalies produced by different-depth density bodies. Results from actual data show that this method extracts the reservoir’s residual density characteristics of plane gravity anomaly on the basis of remaining oil distribution characteristics, consistent with reservoir numerical simulation results. A three-dimensional least-squares inversion of the method for extracting residual gravity anomaly was carried out, with the inversion results consistent with the results of vertical remaining oil distribution characteristics and well-test production results.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Heavy oil is a type of high-viscosity crude oil with high asphaltene and gum content. In China, heavy oil reservoirs are important, widely distributed petroleum resources, among which the Liaohe, Shengli, and Xinjiang oilfields all have large reserves. In order to fully utilize the reservoir and improve oil recovery, it is critical to describe the characteristics of remaining oil in the process of reservoir production and then to formulate a comprehensive reservoir adjustment plan and stimulation measures.

At present, there are relatively few methods to monitor the distribution of remaining oil. Among them, microgravity exploration has become an excellent reservoir-monitoring method, having the advantages of overall monitoring, low cost, no impact on production, nondestructive monitoring, and independence from well verification. With continuous improvement in the accuracy of gravity instruments and continuous progress in technical algorithms, the applications for gravity data are expanding [1]. The gravity exploration method has gradually expanded from understanding regional structural characteristics [2], delineating rock mass range [3], indicating metallogenic prospects [4], seeking local structure [5], and determining stratigraphic rock occurrence [6] to looking for oil and gas resources [7] and describing fluid dynamic changes in oil reservoirs and other fields.

By observing changes in surface gravity data and monitoring changes in oil and gas reservoir density, microgravity exploration can effectively reflect dynamic changes in fluid in the oil and gas reservoir development process and can effectively describe the distribution of remaining oil in heavy oil reservoirs, which is an important link in the dynamic regulation of heavy oil thermal recovery—tapping the potential of remaining oil and prolonging the development life cycle of heavy oil reservoirs. Microgravity-monitoring results for the Alaska Prudhoe Bay gas field have verified that the time-lapse microgravity anomaly can reflect density change caused by water injection in the reservoir and can guide the designs of water injection and enhanced oil recovery [8]. Subsea gravimeters have been used to monitor the height of water-gas contact in water injection gas reservoirs in offshore Norway with an accuracy of meter level [9]. In the Liaohe oilfield, time-shift microgravity monitoring has been used to quantitatively describe changing steam chamber shape during steam-assisted gravity drainage (SAGD) production, providing a reliable basis for adjusting the steam injection scheme [10]. Previous research results have shown that microgravity exploration is suitable for anomalies caused by density changes from displacement and migration of underground fluids in oil and gas reservoirs, with the key solution being to accurately separate the residual gravity anomalies representing reservoir density. At present, the commonly used gravity anomaly separation methods are the filtering method, trend analysis method, peeling method, and nonlinear method [11]. For the filtering method, if the filter is in the nonzero phase, the anomaly extreme point will shift, and there will be significant error when applied to gravity anomaly separation. Affected by the near-source field of the target source, the trend analysis method has difficulty extracting the regional anomaly; the separated local field contains the near-source field information, resulting in anomaly illusion. The premise of the peeling method is to establish a known density model, obtain the anomaly through forward modeling, and then subtract the forward-modeling anomaly from the observed anomaly to obtain the residual anomaly; the accuracy of the established density model directly affects the accuracy of the residual gravity anomaly, and although this method is invalid for regional anomalies, it is effective for near-source effects and high-frequency anomalies. The nonlinear method separates microgravity anomalies and is characterized by multiscale and locality; its disadvantage is asymmetry of surface fitting resulting in distortion of the separated anomaly shape.

In this paper, a multiscale, second-order, surface-fitting, residual gravity anomaly extraction method is proposed. This method combines the multiscale element of the nonlinear separation method and the basic principle of the surface-fitting method. Second-order surface fitting is carried out at different scales to obtain the regional field at the corresponding scale, and then, the calculated regional field is subtracted from the Bouguer gravity anomaly or the regional field at the previous scale to obtain the residual gravity anomaly at this scale. The regional field at different scales can be separated from the residual anomaly and so on, so as to obtain the residual gravity anomaly representing the density of the target area and then explain and analyze the research target. Through forward modeling and field-test verification, the accuracy and reliability of the extraction results of this method were studied in detail.

2. Forward Modeling

Gravity anomaly is the derivative of the additional gravity generated by the residual mass of the geological body to the unit mass at the detection point in the gravity direction. According to the formula of universal gravitation, it can be deduced that the gravity anomaly of a geological body is $$\begin{array}{}\text{(1)}& \Delta g{x}_i,y_i,z_i=G{\iiint }_V\frac{\Delta \rho x_v,y_v,z_v\cdot z_i-z_v}{{{x_i-x_v}^2+{y_i-y_v}^2+{z_i-z_v}^2}^{3/2}}{dx}_v{dy}_v{dz}_v,\end{array}$$where $\Delta g{x}_i,y_i,z_i$ is the derivative of gravity along the $Z$ direction at any coordinate point $x_v,y_v,z_v$, $G$ is the universal gravitation constant, $V$ is the volume of the geological body, and $\Delta \rho x_v,y_v,z_v$ is the residual density difference of a volume element in the geological body coordinate $x_v,y_v,z_v$.

For a reservoir under heavy oil thermal recovery, the buried depth and thickness of the reservoir are known. According to Equation (1), the microgravity value monitored on the surface is mainly controlled by the density change in the reservoir. According to Biot’s theory [12], the original density of reservoir $\rho$ can be calculated by [13] $$\begin{array}{}\text{(2)}& \rho =1-\varphi \cdot 1-V_{\text{sh}}\cdot {\rho }_M+V_{\text{sh}}\cdot {\rho }_{\text{sh}}+\varphi \cdot {\rho }_l,\end{array}$$where $\varphi$ is reservoir porosity, $V_{\text{sh}}$ is the shale content in the reservoir, and ${\rho }_M$, ${\rho }_{\text{sh}}$, and ${\rho }_l$ are rock skeleton density, argillaceous density, and liquid density in the reservoir, respectively.

After heavy oil thermal recovery development, the reservoir density can be expressed by $$\begin{array}{}\text{(3)}& \rho =1-\varphi \cdot 1-V_{\text{sh}}\cdot {\rho }_M+V_{\text{sh}}\cdot {\rho }_{\text{sh}}+\varphi \cdot 1-S_g\cdot {\rho }_l+S_g\cdot {\rho }_g,\end{array}$$where $S_g$ is the saturation after steam injection thermal recovery and ${\rho }_g$ is injected steam density for thermal recovery.

After a period of thermal recovery, the change in reservoir density is expressed by $$\begin{array}{}\text{(4)}& \Delta \rho =\varphi \cdot {\rho }_g-{\rho }_l\cdot S_g.\end{array}$$

The basic principle of microgravity monitoring for a thermal recovery reservoir can be obtained by combining Equations (4) and (1). $$\begin{array}{}\text{(5)}& \Delta g{x}_i,y_i,z_i=G\cdot \varphi \cdot S_g{\iiint }_V\frac{{\rho }_g{x}_v,y_v,z_v-{\rho }_l{x}_v,y_v,z_v\cdot z_i-z_v}{{{x_i-x_v}^2+{y_i-y_v}^2+{z_i-z_v}^2}^{3/2}}{dx}_v{dy}_v{dz}_v.\end{array}$$

Equation (5) shows that for a specific reservoir, the microgravity-monitoring result is directly proportional to the steam saturation in the reservoir porosity. The formula shows the microgravity-monitoring results to be indicative of dynamic changes in fluids in the reservoir.

The Bouguer gravity anomaly obtained from surface observation is the superposition of gravity anomalies generated by all density bodies in the microgravity exploration area (including the gravity anomaly represented by Equation (5)). The gravity anomaly represented by Equation (5) is separated from the Bouguer gravity anomaly by the multiscale, second-order, surface-fitting method. The implementation process is as follows:

Assuming the observed Bouguer gravity anomaly is $gx,y,0$, via the multiscale nonlinear method, its expression can be $$\begin{array}{}\text{(6)}& g_{bg}x,y,0=g_{r1}+g_{l1}=g_{r2}+g_{l1}+g_{l2}=g_{r3}+g_{l1}+g_{l2}+g_{l3}\cdots ,\text{(7)}& {\Delta g_r}_i=g_{bg}-{g_r}_i=\sum _{i=1}^n{g}_{li},\end{array}$$where $x,y,0$ represent the ground coordinates, $g_{l1},g_{l2},g_{l3},\cdots$ represent local anomalies of different scales, and ${\Delta g_r}_i$ is the residual gravity anomaly at the $i$^th scale, where $i=1,2,3,4$.

The local fields ${g_r}_i$ at different scales are obtained by the second-order surface-fitting method. Located in a small area centered on a point $m,n,0$, $g_r$ can be expressed as $$\begin{array}{}\text{(8)}& g_r=a_0+a_{1}x+a_{2}y+a_{3}xy+a_4{x}^2+a_5{y}^2,\end{array}$$where $x,y$ is the distance from the point in the area to the center point. Then, the sum of error squares at each point is expressed as $$\begin{array}{}\text{(9)}& Eg,g_r=\sum _i\sum _j{gi,j-g_{r}i,j}^2,\end{array}$$where $i$ and $j$ are taken in the area with the point as the center $m,n$.

Firstly, the second-order surface coefficient $a_0,a_1,a_2,a_3,a_4,a_5$ is obtained by the least-squares method, and then, the same operation is carried out for each point to obtain the regional field value of each point under the corresponding scale. Next, the residual gravity anomaly under the corresponding scale is calculated by Equation (7), and finally, the reasonable scale is determined according to the detection target depth and relevant geological data. The residual gravity anomaly generated by density body in the depth range of microgravity detection target is obtained.

2.1. Modeling

In order to verify the method, a model for identifying the effect of gravity separation was adopted in this study. The model was first proposed by Guo et al. [14]. The model parameters are shown in Table 1, and the model distribution is shown in Figure 1.

Table 1

Forward-modeling parameters.

No.	Length (m)	Width (m)	Thickness (m)	Depth of top (m)	Residual density (g/cm3)
A1	4000	5000	1500	1500	0.15
A2	3000	7000	1500	1500	−0.10
A3	1500	3300	1500	1500	0.10
B1	400	1600	400	400	−0.25
B2	400	1000	400	400	0.25
B3	600	800	400	400	0.25
B4	600	900	400	400	0.25
B5	600	2400	400	400	−0.25
C1	100	200	100	100	0.50
C2	100	200	100	100	−0.50
C3	100	100	100	100	0.50

[figure(s) omitted; refer to PDF]

The model is divided into 3 classes (A, B, and C) composed of 11 cuboids at different depths. The abnormal body with positive residual density simulates the deposits with high density such as metal ore, while the negative residual density simulates the karst cave, the steam cavity, and the abnormal body with low density and high porosity. For the abnormal bodies, class A represents deep background fields, such as basement fluctuation and background density body deeper than the target layer; class B represents the studied target body, which is an effective signal to be retained in the separation process; class C is removed as a shallow interference signal in some cases and is retained as a research target in some cases, such as shallow steam-channeling research.

For the model, an observation system with a sampling interval of 10 m and survey network of $2001\times 2001$ is used for forward-modeling simulation. The forward-modeling results are shown in Figure 2(a).

[figure(s) omitted; refer to PDF]

For simulating the distribution of gravity fields, random Gaussian noise with a standard deviation of $0.2\times {10}^{-5}$ m/s² is added to the theoretical gravity anomaly (Figure 2(b)). Figure 2 shows the superposition results of anomalies generated by various geological bodies at different depths, testing the feasibility of the anomaly extraction method in depicting, distinguishing, and multiscaling local details.

2.2. Separated Result Analysis of Forward Modeling

The multiscale, second-order, surface-fitting method proposed in this paper extracts the residual gravity anomalies of the gravity field generated by the model at different scales (2ⁿ times the sampling interval). The calculation results according to Equation (7) are shown in Figure 3: among the residual gravity anomalies separated at different scales, the anomalies generated by various geological bodies are gradually displayed with the continuous increase of scales; they are consistent with geological body location.

[figure(s) omitted; refer to PDF]

In previous studies, Guo et al. [14] and Shi et al. [15] used a variety of methods to carry out anomaly separation experiments on the forward gravity field of the model. It was found that although these methods can extract the residual gravity anomalies generated by the model geological body, they cannot distinguish the residual gravity anomalies generated by the geological body at different depths—and not at multiscale. The multiscale surface-fitting method can effectively separate the regional field at different depths from the residual field without distortion of local anomalies. The anomaly extraction results are more accurate and amplitude-preserved than other methods.

From the analysis of forward-modeling results, it can be seen that the multiscale, second-order, surface-fitting gravity anomaly extraction method proposed in this paper has the following advantages:

(1) It is multiscale and supplies local description of anomalies without introducing abnormal distortion along the $X$ and $Y$ directions.

(2) The residual gravity anomaly of each scale can reflect the gravity effect produced by the abnormal body at different depth levels. The scale can be selected according to the exploration target, and the abnormal position is clear and accurate.

(3) It has a strong ability to suppress random noise without involving any empirical formula or parameters in anomaly extraction. Therefore, it is not susceptible to subjective influence and has a good theoretical basis and experimental effect.

3. Field Application

In 2019, microgravity monitoring was carried out in a heavy oil reservoir development area in Xinjiang. Using the multiscale, second-order, surface-fitting method proposed in this paper, the residual gravity anomaly representing the density body in the reservoir was extracted according to the scale determined by the buried depth range of the exploration target, and the residual oil distribution was described in slices. In addition, based on the gravity anomaly generated by the reservoir density volume extracted by multiscale surface fitting, the reservoir density volume was obtained using the three-dimensional (3D) least-squares inversion method, and the remaining oil distribution was characterized vertically. Finally, the horizontal and vertical characterization results of remaining oil distribution were verified by well data and production performance data.

3.1. Horizontal Characterization of Remaining Oil Distribution

The residual gravity anomaly results obtained by the multiscale, second-order, surface-fitting method representing the density in the reservoir are shown in Figure 4. On the whole, the warm-color area with a relatively high amplitude of residual gravity anomaly indicates relatively high reservoir density, high remaining oil saturation, and significant remaining oil development potential. The cold-color area with a relatively low amplitude of residual gravity anomaly indicates relatively low reservoir density, low steam sweep, low saturation of remaining oil, and low remaining oil development potential.

[figure(s) omitted; refer to PDF]

In order to verify the depiction results of microgravity monitoring on the plane, the reservoir density was simulated and calculated using the numerical simulation method, combined with the actual production performance data, logging, and other basic data, after which the weighted average result of reservoir density by depth was obtained. In Figure 5, the warm-color area with a relatively high density value represents high residual oil saturation. On the contrary, the cold-color area with a relatively low-density value represents low residual oil saturation. Comparing Figures 4 and 5, the distribution results from microgravity monitoring for the areas with high and low residual oil saturations are consistent with the residual oil distribution results obtained by the numerical simulation method, verifying the residual oil distribution results from microgravity monitoring on the plane.

[figure(s) omitted; refer to PDF]

Therefore, the distribution of residual oil above different sections of the horizontal well can be described according to the microgravity-monitoring results, as shown in Figure 4, where the red section of the horizontal well indicates the reservoir above the well section having high residual oil saturation, the black section indicates medium residual oil saturation, and the white section indicates low residual oil saturation.

3.2. Vertical Characterization of Remaining Oil Distribution

Figure 6 shows the density body section along the horizontal well direction. On the whole, the reservoir density is shown to gradually decrease from the horizontal well position to the top of the reservoir, indicating a high recovery degree of the upper reservoir and the remaining oil being mainly distributed in the middle and lower parts of the reservoir. Along wells HW22, HW23, and HW24, the remaining oil is shown to be mainly distributed in the middle section, and along well HW25, the remaining oil is shown to be mainly distributed in the middle and end sections.

[figure(s) omitted; refer to PDF]

Figure 7 shows density body sections along the source direction (Figure 7(a)) and perpendicular to the source direction (Figure 7(b)). The comparison shows the reservoir density between wells along the material source direction being relatively low, indicating the degree of remaining oil production as high and the connectivity along the material source direction being better.

[figure(s) omitted; refer to PDF]

As can be seen from Figures 6 and 7, the density of the reservoir decreases from bottom to top. The remaining oil is mainly distributed in the middle and lower parts of the reservoir, which are macroscopically antirhythmic in the vertical direction. At the same time, the remaining oil is affected by steam overlap, resulting in quality deficit and density reduction in the upper part of the reservoir, highly corresponding to production injection ratio and low remaining oil content. The steam sweep degree in the lower part of the reservoir is low, and oil and water accumulate in the lower part of the reservoir under the effect of gravity, which is characterized by increased mass, increased density, low production injection ratio, and high residual oil.

Well data verified the characterization results of residual oil distribution in the vertical direction from microgravity monitoring. First, temperature tests on well 192 and well 234 were used for well temperature monitoring. The peak temperature area is the steam chamber development area; the reservoir density in this area should be low. The inversion density body sections are along well 192 and well 234 (Figure 8). It can be seen that the depth range of the low-density area of the inversion density body along the well is consistent with the depth range of the low-density area of the corresponding well temperature monitoring—verifying the accuracy of extracting and characterizing the residual gravity anomaly generated by reservoir density using the multiscale, second-order, surface-fitting method.

[figure(s) omitted; refer to PDF]

Then, the production performance data were used to further verify the description results of residual oil distribution in the vertical direction from microgravity monitoring. The statistics of all production performance data from the microgravity exploration block show the recovery degrees of J₃q₂^2-1 and J₃q₂^2-2 layers to be 26.4% and 19.7%, respectively, indicating low remaining oil saturation in the upper part of the reservoir and high remaining oil saturation in the middle and lower parts. The microgravity-monitoring results for remaining oil distribution in the vertical direction are consistent with the production performance results. The single-well production performance statistics show the reservoir area with the low-production-injection-ratio well to have a relatively high reservoir density and relatively high residual oil saturation. On the contrary, the reservoir area with the high-production-injection-ratio well is shown to have a relatively low corresponding reservoir density and relatively low remaining oil saturation. It can be seen that the microgravity-monitoring results are consistent with the statistical results of single-well production performance data (Figure 9)—further verifying the accuracy of the residual gravity anomaly characterized by reservoir density extracted by multiscale, second-order surface fitting.

[figure(s) omitted; refer to PDF]

4. Conclusions

(1) In this paper, a multiscale, second-order, surface-fitting separation method of gravity anomalies is proposed; it is suitable for separating microgravity field anomalies. This nonlinear method is more sensitive to abnormal signal changes and does not cause abnormal distortion. The anomaly extraction result is accurate, and the amplitude remains unchanged.

(2) Compared with the traditional gravity anomaly separation method, the proposed method can effectively separate the gravity anomalies of anomaly bodies at different depths and can accurately describe the spatial distribution and edge characteristics of anomaly bodies.

(3) The separation method was applied for a field test of microgravity residual oil characterization. The gravity anomaly representing the residual density of the reservoir extracted based on the multiscale, second-order, surface-fitting method can effectively characterize the distribution of remaining oil in the horizontal and vertical directions. The reservoir numerical simulation results and production performance data verify the reliability of the results, and the 3D density distribution provides a reliable basis for the development of remaining oil.

Acknowledgments

This work was supported by the China National Petroleum Corporation Limited Upstream Field Forward-looking Basic Technology Research Project.

References

[1] A. L. Codd, L. Gross, A. Aitken, "Fast multi-resolution 3D inversion of potential fields with application to high-resolution gravity and magnetic anomaly data from the Eastern Goldfields in Western Australia," Computers & Geosciences, vol. 157, article 104941,DOI: 10.1016/j.cageo.2021.104941, 2021.

[2] A. V. Satyakumar, A. K. Pandey, A. P. Singh, M. Tiwari Virendra, "Delineation of structural and tectonic features in the Mahanadi basin, eastern India: new insights from remote sensing and land gravity data," Journal of Asian Earth Sciences, vol. 227, article 105116,DOI: 10.1016/j.jseaes.2022.105116, 2022.

[3] M. B. Dejene, F. Yasuhiro, N. Jun, S. Hakim, "Interpretation of gravity data to delineate the subsurface structures and reservoir geometry of the Aluto-Langano geothermal field, Ethiopia," Geothermics, vol. 94, article 102093,DOI: 10.1016/j.geothermics.2021.102093, 2021.

[4] K. Hassan, A. Firouz, G. Abdolreza, "Regional magnetic and gravity structures and distribution of mineral deposits in Central Iran: implications for mineral exploration," Journal of Asian Earth Sciences, vol. 217, article 104828,DOI: 10.1016/j.jseaes.2021.104828, 2021.

[5] B. U. Camilus, A. A. Effiong, E. E. Ohara, A. Emmanuel, "Novel technique for the interpretation of gravity anomalies over geologic structures with idealized geometries using the Manta ray foraging optimization," Journal of Asian Earth Sciences: X, vol. 6, article 100070,DOI: 10.1016/j.jaesx.2021.100070, 2021.

[6] Y. Abdelfettah, J.-J. Tiercelin, P. Tarits, S. Hautot, M. Maia, P. Thuo, "Subsurface structure and stratigraphy of the northwest end of the Turkana Basin, Northern Kenya Rift, as revealed by magnetotellurics and gravity joint inversion," Journal of African Earth Sciences, vol. 119, pp. 120-138, DOI: 10.1016/j.jafrearsci.2016.03.008, 2016.

[7] H. Aghajani, A. Moradzadeh, H. Zeng, "Detection of high-potential oil and gas fields using normalized full gradient of gravity anomalies: a case study in the Tabas Basin, Eastern Iran," Pure and Applied Geophysics, vol. 168 no. 10, pp. 1851-1863, DOI: 10.1007/s00024-010-0169-y, 2011.

[8] J. L. Hare, J. F. Ferguson, J. L. Brady, "The 4D microgravity method for waterflood surveillance: part IV — modeling and interpretation of early epoch 4D gravity surveys at Prudhoe Bay, Alaska," Geophysics, vol. 73 no. 6, pp. WA173-WA180, DOI: 10.1190/1.2991120, 2008.

[9] M. Zumberge, H. Alnes, O. Eiken, G. Sasagawa, T. Stenvold, "Precision of seafloor gravity and pressure measurements for reservoir monitoring," Geophysics, vol. 73 no. 6, pp. WA133-WA141, DOI: 10.1190/1.2976777, 2008.

[10] Y. Li, F. Ren, L. Yang, D. Zhou, X. Tian, "Description of steam chamber shape in heavy oil recovery using 4D microgravity measurement technology," Petroleum Exploration and Development, vol. 40 no. 3, pp. 409-412, DOI: 10.1016/S1876-3804(13)60052-X, 2013.

[11] O. A. Stepanov, D. A. Koshaev, A. V. Motorin, "Identification of gravity anomaly model parameters in airborne gravimetry problems using nonlinear filtering methods," Gyroscopy and Navigation, vol. 6 no. 4, pp. 318-323, DOI: 10.1134/S2075108715040136, 2015.

[12] M. A. Biot, "Theory of propagation of elastic waves in a fluid-saturated porous solid," The Journal of the Acoustical Society of America, vol. 28 no. 2, pp. 168-178, DOI: 10.1121/1.1908239, 1956.

[13] L. Yong, L. Zhengwen, L. Zhirong, L. Qiong, "Study on reservoir density prediction technology," Oil Geophysical Prospecting, vol. 42 no. 2, pp. 216-219, 2007.

[14] L. H. Guo, X. H. Meng, L. Shi, Z. X. Chen, "Preferential filtering method and its application to Bouguer gravity anomaly of Chinese continent," Chinese Journal of Geophysics, vol. 55 no. 12, pp. 4078-4088, 2012.

[15] S. Lei, C. Shi, J. C. X. Weimin, H. Lu, F. Guo, "Analysis the of characteristics of gravity anomaly in the seismic area of Lushan earthquake based on preferential filtering method," Earthquake Science, vol. 35 no. 5, pp. 704-716, 2013.