Numeric Simulation of Acoustic-Logging of Cave Formations

Abstract

The finite difference (FD) method of monopole source is used to simulate the response of full-wave acoustic-logging in cave formations. The effect of the cave in the formation of borehole full-waves was studied. The results show that the radius of cave is not only linearly related to the first arrival of the compressional wave (P-wave), but also to the energy of the shear wave (S-wave). The converted S (S–S wave) and P-waves (S–P wave) are formed when the S-wave encounters the cave. If the source distance is small, the S–S and S–P waves are not separated, and the attenuation of the S-wave is not large, due to superposition of the converted waves. The S–P wave has been separated from the S-wave when the source distance is large, so the attenuation of the S-wave increases. The amplitude of the P and S–waves changes most when the distance of the cave to the borehole wall reaches a certain value; this value is related to the excitation frequency. The amplitude of the Stoneley wave (ST wave) varies directly with the radius of cave. If the radius of the cave is large, the energy of ST wave is weak. The scattered wave is determined by the radius and position of the cave. The investigation depth of a monopole source is limited. When the distance of the cave to the borehole wall exceeds the maximum investigation depth, the borehole acoustic wave is little affected by the cave. In actual logging, the development of the cave can be evaluated by using the first arrival of the P-wave and the energy of the S and ST waves.

1. Introduction

Caves are widely distributed in carbonate and igneous reservoirs, causing heterogeneity and making it difficult to evaluate the reservoirs by conventional logging. Based on the observation of cores, the reservoir spaces of igneous and carbonate formations can be divided into primary pores and secondary pores, according to their genesis [1]. Secondary pores contain caves. These can be divided into matrix-dissolution pores, intragranular-dissolution pores and apricot-kernel dissolution pores in igneous formations [2]. According to the genesis, carbonate formations have porous caves, intergranular caves and fractured caves. On the basis of the size of the cave, they may fall into middle and small-type caves and large caves. Generally, caves with a diameter greater than 10 mm can be regarded as large caves [3]. Because dissolution pores and caves often transform primary pores, making them better connected and improving the reservoir property of the rock, they are important reservoir spaces [4]. In fact, the caves in the formation can be defined as voids, washouts or large pores of secondary porosity in a more general. The existence of caves greatly changes the speed and energy of acoustic-logging [5]. Therefore, it is very important to grasp the effect of caves on borehole acoustic wave propagation, which plays an important role in the detection and evaluation of caves. Currently, the effect of fractured and gas-bearing reservoirs on borehole acoustics has been investigated [6,7,8]. However, there are few studies on the acoustic-logging response of caves. Heterogeneous bodies such as fractures and caves, scatter seismic-wave propagation [9]. Elastic-wave scattering in heterogeneous media can be divided into “velocity” and “wave impedance” types [10]. Due to the fact that the scattered wave of a borehole dipole source is mainly a SH shear wave [11,12,13], the dipole source has been used to study the scattering effect in heterogeneous formation [14]. The relationship between S-wave wavelength and the size of cave can be obtained by dipole acoustic-logging [15]. According to the different filling materials, caves can be divided into four types: unfilled, sand-mud-filled, breccia filled, and calcite-filled [16]. The filling degree of a cave can slow P and S-waves. Studies have shown that different filling degrees of caves has a good corresponding relationship with acoustic simulation curves, and the results have certain application effect [17]. In cave formation, the slowness of the P-wave increases and the density of well-logging decreases. The resistivity is generally low in conventional logging [18]. Acoustic-logging and electrical-imaging logging can evaluate caves comprehensively; ST wave energy of acoustic-logging reduces due to caves [3].

The above studies show that it is feasible to use acoustic-logging to evaluate near-borehole caves. Based on actual logging data, the results have qualitatively shown the effects of cave on conventional logging, such as P-wave slowness. The response of component waves to caves has been studied. There is no detailed regularity of the size and position of cave on borehole acoustic waves. Monopole is a common source in acoustic-logging [19,20], and the component waves with different frequencies formed by monopole contain the information of specific formations. It is difficult to analyze scattered-wave fields of monopole acoustic-logging because of their complexity. Numeric simulation in cave reservoirs by monopole sources has not been reported. In this study, the finite difference (FD) method of the elastic wave equation is used to simulate a monopole full-wave acoustic-logging response of a fluid-filled cave in a homogeneous formation [21,22]. On this basis, the results of a numeric simulation are used to analyze actual logging data [23].

In the process of simulating a wave field inside and outside of a borehole by FD, matrix calculations are mainly carried out. It is necessary to divide small grids to meet the calculation requirements of cave, which greatly increases the time costs of the numeric simulation. Therefore, computer-unified device architecture (CUDA) parallel-computing was used to improve the calculation efficiency [24,25,26,27]. The essence of CUDA is to change the computing mode from central processing unit (CPU) serial to graphics-processing unit (GPU) parallel.

2. Theory and Methods

2.1. FD Method for Elastic Wave Equation

In order to study the issue more simply, we use the elastic wave equation in the 2D-cylindrical coordinate system. A 2D-model simulation has advantages over a 3D-model simulation in terms of calculation times. A 2D calculation of a cave can also show the main characteristics of the problem studied [15]. The staggered-grid finite difference method is used to simulate the propagation of elastic waves in our simulation [6]. The wave equation of velocity stress in isotropic media in cylindrical coordinates is as follows:

$$\begin{gathered} \left[{\tau }_{rr\left(i+\frac{1}{2},j+\frac{1}{2}\right)}^{n+1}-{\tau }_{rr\left(i+\frac{1}{2},j+\frac{1}{2}\right)}^n\right]/\mathsf{\Delta }t=\left[{\lambda }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}+2{\mu }_{i+\frac{1}{2},j+\frac{1}{2}}\right]\left[L_R^1{v}_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}\right] \\ +{\lambda }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}\left[L_Z^1{v}_{z\left(i+\frac{1}{2},j\right)}^{n+\frac{1}{2}}+{\sigma }_i{v}_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}/r_{i+\frac{1}{2}}\right] \end{gathered}$$

(1)

$$\begin{gathered} \left[{\tau }_{zz\left(i+\frac{1}{2},j+\frac{1}{2}\right)}^{n+1}-{\tau }_{zz\left(i+\frac{1}{2},j+\frac{1}{2}\right)}^n\right]/\mathsf{\Delta }t={\lambda }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}\left[L_R^1{v}_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}+\frac{{\sigma }_i{v}_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}}{r_{i+\frac{1}{2}}}\right] \\ +\left[{\lambda }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}+2{\mu }_{i+\frac{1}{2},j+\frac{1}{2}}\right]\left[L_Z^1{v}_{z\left(i+\frac{1}{2},j\right)}^{n+\frac{1}{2}}\right] \end{gathered}$$

(2)

$$\begin{gathered} \left[{\tau }_{\theta \theta \left(i+\frac{1}{2},j+\frac{1}{2}\right)}^{n+1}-{\tau }_{\theta \theta \left(i+\frac{1}{2},j+\frac{1}{2}\right)}^n\right]/\mathsf{\Delta }t={\lambda }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}\left[L_R^1{v}_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}+L_Z^1{v}_{z\left(i+\frac{1}{2},j\right)}^{n+\frac{1}{2}}\right] \\ +\left[{\lambda }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}+2{\mu }_{i+\frac{1}{2},j+\frac{1}{2}}\right]\left[{\sigma }_i{v}_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}/r_{i+\frac{1}{2}}\right] \end{gathered}$$

(3)

$$\left[{\tau }_{rz\left(i,j\right)}^{n+1}-{\tau }_{rz\left(i,j\right)}^n\right]/\mathsf{\Delta }t={\mu }_{\left(i-\frac{1}{2},j-\frac{1}{2}\right)}^H\left[L_R^1{v}_{z\left(i-\frac{1}{2},j\right)}^{n+\frac{1}{2}}+L_Z^1{v}_{r\left(i,j-\frac{1}{2}\right)}^{n+\frac{1}{2}}\right]$$

(4)

$$\begin{array}{cc}\hfill {\sigma }_i{\rho }_{\left(i-\frac{1}{2},j+\frac{1}{2}\right)}\left[v_{r\left(i,j+\frac{1}{2}\right)}^{n+\frac{1}{2}}-v_{r\left(i,j+\frac{1}{2}\right)}^{n-\frac{1}{2}}\right]/\mathsf{\Delta }t=& L_R^1{\tau }_{rr\left(i-\frac{1}{2},j+\frac{1}{2}\right)}^n+L_Z^1{\tau }_{rz\left(i,j\right)}^n\hfill \\ & +\left[{\sigma }_i{\tau }_{rr\left(i-\frac{1}{2},j+\frac{1}{2}\right)}^n-{\sigma }_i{\tau }_{\theta \theta \left(i-\frac{1}{2},j+\frac{1}{2}\right)}^n\right]/r_i\hfill \end{array}$$

(5)

$${\sigma }_i{\rho }_{\left(i+\frac{1}{2},j-\frac{1}{2}\right)}\left[v_{z\left(i+\frac{1}{2},j\right)}^{n+\frac{1}{2}}-v_{z\left(i+\frac{1}{2},j\right)}^{n-\frac{1}{2}}\right]/\mathsf{\Delta }t=L_R^1{\tau }_{xz\left(i,j\right)}^n+L_Z^1{\tau }_{zz\left(i+\frac{1}{2},j-\frac{1}{2}\right)}^n+{\sigma }_i{\tau }_{rz\left(i,j\right)}^n/r_{i+\frac{1}{2}}$$

(6)

In the equations mentioned above, the subscripts i and j in brackets refer to spatial index of nodes in x and z directions, respectively, the superscript n represents the time discrete index, ${\tau }_{rr}$, ${\tau }_{\theta \theta }$ and ${\tau }_{zz}$ refer to normal stress in x, $\theta$ and z directions, ${\tau }_{rz}$ is shear stress, $\mathsf{\Delta }t$ is the time step, $L_R^1$ and $L_Z^1$ are second-order difference operators of space, and the parameters ρ and μ represent the density and shear moduli of the medium. On the medium interface, the arithmetic mean value of ρ and the harmonic mean value of μ are employed; on the solid–liquid interface, μ is set to 0 [6,7]. ${\sigma }_i{f}_{\left(i,j\right)}$ and ${\sigma }_j{f}_{\left(i,j\right)}$ are arithmetic average operators and defined as:

$${\sigma }_i{f}_{\left(i,j\right)}=0.5\left[f_{\left(i+1,j\right)}+f_{\left(i,j\right)}\right]$$

(7)

$${\sigma }_j{f}_{\left(i,j\right)}=0.5\left[f_{\left(i,j\right)}+f_{\left(i,j+1\right)}\right]$$

(8)

${\mu }_{\left(i,j\right)}^H$ is the harmonic mean and defined as

$${\mu }_{\left(i,j\right)}^H=4/\left[1/{\mu }_{\left(i+\frac{1}{2},j+\frac{1}{2}\right)}+1/{\mu }_{\left(i+\frac{1}{2},j-\frac{1}{2}\right)}+1/{\mu }_{\left(i-\frac{1}{2},j+\frac{1}{2}\right)}+1/{\mu }_{\left(i-\frac{1}{2},j-\frac{1}{2}\right)}\right]$$

(9)

To reduce the numeric dispersion, the space grid needs to satisfy the following condition:

$$\mathsf{\Delta }\mathrm{r}\le v_{min}/10f_{max}$$

(10)

where $v_{min}$ is the minimum velocity of numeric model, and $f_{max}$ is the maximum frequency of the source [6,7]. The stability condition of the second-order FD is

$$v_{max}\mathsf{\Delta }t{\left[\frac{1}{\mathsf{\Delta }{r}^2}+\frac{1}{\mathsf{\Delta }{z}^2}\right]}^{\frac{1}{2}}<1$$

(11)

where $v_{max}$ is the maximum velocity of the numeric model [6,7]. The monopole source function is

$$s\left(t\right)=\{\begin{array}{c}\frac{1}{2}\left[1+\cos \frac{2\pi }{T_c}\left(t-\frac{T_c}{2}\right)\right]cos2\pi f_0\left(t-\frac{T_0}{2}\right)，0\le t\le T_c\\ 0，t>T_c\end{array}$$

(12)

where $f_0$ denotes the center frequency of the source, t represents the time, and $T_c$ is the acoustic source pulse width, $T_c=2/f_0$.

2.2. The Principle of CUDA Parallel-Computing

The CUDA parallel-computing model of FD method of elastic wave equation is established; the algorithm implementation flow is shown in Figure 1. During the initialization of the calculation model, the CPU first allocates storage space in the GPU memory for the velocities, stresses and medium parameters. The kernel functions of stress, velocity fields and source are then designed. The GPU is controlled by the CPU to complete the step-by-step calculations of velocity and stress fields. It was found that for the homogeneous medium, the calculation time of the conventional program was 10,571.5 s; the calculation time of the parallel program was only 704.7 s through the simulation experiment. CUDA parallel-computing is used as a computing method in this study, so it is not introduced in detail here.

2.3. Borehole Model and Non-Splitting Perfect-Matched Layer Absorbing Boundary Condition

Figure 2 shows a borehole model surrounded by an elastic formation. The cave was filled with water and the formation was homogeneous. Perfect-matched layers (PML) were added at the top, bottom and right side of the model on the $\mathrm{r}\_\mathrm{z}$ plane. The model dimension was 4 m × 2 m, the grid step was $\mathsf{\Delta }r=\mathsf{\Delta }z=2\mathrm{mm}$, and the time step was $\mathsf{\Delta }t=0.2\mathsf{\mu }\mathrm{s}$. A non-splitting perfect-matched layer (NPML) approach [28] was used to mitigate the ghost reflection at the artificial boundary of the model.

The elastic parameters of each medium are shown in

Table 1. Physical parameters of different media.

Elastic Parameters	${\mathit{v}}_{\mathit{p}}\mathbf{m}\mathbf{/}\mathbf{s}$	${\mathit{v}}_{\mathit{s}}\mathbf{m}\mathbf{/}\mathbf{s}$	$\mathit{\rho }\mathbf{k}\mathbf{g}\mathbf{/}{\mathbf{m}}^{\mathbf{3}}$
Elastic formation	3670	2170	2650
Borehole fluid	1500	0	1000
Cave fluid	1500	0	1000

. Alteration usually occurs in carbonate and igneous reservoirs. The rock density of this kind of reservoir is generally between 2.5–2.8 $\mathrm{g}/{\mathrm{cm}}^3$. In this study, 2.65 $\mathrm{g}/{\mathrm{cm}}^3$ was selected as the density of the formation [20]. The P and S-wave velocity were selected according to the principle that the P-wave velocity is generally 1.6–1.9 times that of the S-wave velocity [20]. The selection of these physical parameters is similar in numeric simulations of fractures and cave formations [7,15,29]. If the physical parameters of calcite and dolomite are selected in the numeric simulation, it will be of more practical significance.

In order to verify the correctness of the results of the FD method, the propagation of the ST wave in homogeneous elastic formation was simulated by using the FD and the RAI method [29]. The formation parameters are shown in Table 1. The center frequency of the source was 2 kHz. Figure 3 shows the comparison between the waveform calculated by the FD and RAI method. The waveforms of the two methods matched well and it showed that the FD method was effective. In addition, there was a little difference in waveform amplitude between the two methods, which was caused by discrete Fourier transform (DFT) in RAI method and finite difference dispersion in FD method. The maximum error of the two methods was 1.2%.

3. Numeric Simulation Results

3.1. Borehole Acoustic Wave and Wave Field Snapshot of Cave Formation

Figure 4 shows the full-wave waveform of cave formation. In this study, a monopole source is selected. The P and S-wave cannot be excited when the center frequency of the source is 2 kHz. There is only a ST wave in the borehole. The ST wave is unchanged because the cave in the formation cannot directly affect it in Figure 4a. The P, S and pseudo-Rayleigh waves appear in the wave train (Figure 4b–d) with the increase of the source frequency. Scattering waves are formed when the P and S-waves encounter the water-filling cave in the formation. Borehole acoustic waves are also changed because of the superposition with scattered waves [30]. In order to suppress pseudo-Rayleigh wave, the center frequency of source is usually 8–12 kHz.

The center frequency of monopole source was 10 kHz in this study. Figure 5. shows the wave field of acoustic wave at different times in cave formation. The acoustic wave propagated along the borehole or formation to the PML boundary of the model. There was no artificial boundary reflection, and this proves that the NPML was effective. The P-wave encountered the cave at 0.9 ms. Because the scattered wave formed by P-wave was weak, only the original wave field excited by source can be seen in Figure 5a. The S–P and S–S waves were formed when an S-wave encountered a cave due to energy conversion. The S–P wave can be observed in Figure 5b beside the original wave field formed by the source at 1.3 ms. The scattered-wave field formed by the S–wave can be clearly observed in Figure 5c—including S–P and S–S waves. It should be noted that the amplitude of the S-wave decreased significantly after it passed through the cave in Figure 5c (the color of the S-wave became lighter).

3.2. The Effect of the Cave on the Borehole P-Wave

Figure 6a shows the effect of different radii of cave on P-waves.

• The amplitude of the P-wave varied slightly with the radius of the cave;
• The larger the cave radius was, the greater was the delay of the first arrival of the P-wave (the first arrival in our study was defined as the maximum of the first phase). From Figure 7a, the relationship between the first arrival of the P-wave and the radius of the cave was linear;
• Figure 6b shows the effect of distance of the cave to the borehole wall on the P-wave;
• There exists a distance of the cave to the borehole wall defined as $L_P$. When the distance was equal to $L_P$, the amplitude of the P-wave generally decreased in observations. Here $L_P$ = 16 cm in our model;
• It was found that $L_P$ was related to the excitation frequency according to the test. $L_P$ = 24 cm if the excitation frequency was 7 kHz in Figure 7c and $L_P$ = 16 cm if the excitation frequency was 10 kHz in Figure 7b;
• When the distance was more than 24 cm (the excitation frequency was 10 kHz), it had little effect on the P-wave because of the limitation of investigation depth of the monopole source in our model.

Because the cave was fully filled with water and the wave propagation speed in the fluid was less than that in the formation, the first arrival of the P-wave was delayed. The larger the radius of the cave, the longer was the P-wave propagation time in the fluid, and the later the first arrival time of the P-wave. In this study, a single cave was simulated. There may be multiple caves of different sizes in an actual formation, which will cause more delays on the first arrival of the P-wave [31]. P-wave recorded by the receiver was the superposition of scattered wave and “gliding” P-wave refracted on the borehole wall. The scattered wave formed by the cave changed with different distances to the borehole wall (L), so the amplitude of the P-wave may be different when the scattered wave is superimposed with it. In this study, the P-wave amplitude was reduced slightly, and it may increase if the source distance was changed because the scattered wave also varied with different source distance.

3.3. The Effect of the Cave on Borehole S-Wave

Figure 8a shows the effect of the radius of the cave on the S-wave.

• The attenuation of the S-wave was large (the first arrival in our study was defined as the maximum of the first phase). With the increase of cave radius, the S-wave attenuation became stronger;
• The energy of the S-wave was more sensitive to cave radius than P-wave;
• Figure 8b shows the effect of the distance of the cave on S-wave;
• There exist a distance of cave to the borehole wall defined as $L_S$. When the distance was equal to $L_S$, the amplitude of the S-wave decays most by observation. Here $L_S$ = 16 cm in our model. $L_S$ was also related to the excitation frequency;
• When the distance was more than 24 cm (the center frequency of source was 10 kHz), it had little effect on the S-wave in our model. Because of the investigation depth of the monopole source was limited. A distance greater than 24 cm may exceed the maximum investigation depth.
• S-wave was more sensitive to the change of the distance of the cave than P-wave.

When the S-wave encounters the water-filled cave, the water inhibits the propagation of the S-wave and a considerable part of the energy of the S-wave is converted into the compressed energy (the energy of S–P wave) of the fluid [31]. The compression energy may be converted into S-wave and other mode waves, and these converted waves are recorded by the receiver finally. Therefore, the S-wave energy recorded by the receiver above the cave will be partially attenuated due to scattering and conversion.

3.4. The Relationship between the Attenuation Coefficient of S-Wave and the Cave

It is found that the existence of a cave can make the S-wave energy decay through simulation and analysis. In this study, the attenuation coefficient A of the S-wave under the action of the cave is defined as

$$A=1-\frac{A_2}{A_1}$$

(13)

where $A_1$ and $A_2$ are the amplitude of the S-wave for a homogeneous formation and a cave formation in the same source distance, respectively. The dotted line is a fitting linear equation in Figure 9a, and there is a positive linear correlation between the attenuation coefficient and the radius of the cave. The slope of linear equation of attenuation coefficient is variable when the distance of cave is different. The large slope indicates that the attenuation coefficient is sensitive to the radius of the cave. When L = 4 cm, the attenuation coefficient is the most sensitive to the cave radius in our model.

When the S-wave propagates in the formation and run into the water-filled cave, it will form S–P and S–S waves. The S–P and S–S waves recorded by the receiver near the cave are superposed together, so the attenuation coefficient is smaller as shown in Figure 9b. The S–P and S–S waves are separated gradually with the increase of the source distance. When the S–P leaves the S-wave, it takes its energy away from the S-wave, so the attenuation coefficient increases [31]. The effect of the cave on the S-wave is reduced with the increase of the source distance, so the attenuation coefficient decreases gradually. In our model, the S–S wave was separated completely when the source distance is 2.6 m as can be seen from Figure 9b.

3.5. The Effect the of Cave on Stoneley Wave

The P and S-wave which can propagate in the formation are formed when the center frequency of source is high. At this time, the scattered wave propagates to the borehole and affects the ST wave. Figure 10a is the wave train in homogeneous elastic formation. Figure 10b–f represent the wave train in the formation with cave of different radii. It can be seen that when the radius of the cave is small (Figure 10b,c), the attenuation of the ST wave is weak compared with Figure 10a. When the radius of the cave is large (Figure 10d–f), the amplitude of the ST wave has a great attenuation. Therefore, the attenuation of the ST wave not only corresponds to fractures, but also to caves.

Figure 11a shows acoustic waveforms. In order to study the relationship between the ST wave and the cave radius, it is necessary to calculate the energy of the ST wave. The full-wave frequency spectrum is obtained by Fourier transform in Figure 11b; the low-pass filter is used to obtain ST wave. The energy of the ST wave is obtained by summing the square of the corresponding value of each frequency. The frequency of the ST wave is less than 4 kHz according to Figure 11b. The cutoff frequency of the low-pass filter is 4 kHz. The black line in Figure 11a is the low-pass filtered waveform. The S-wave and pseudo-Rayleigh wave have been suppressed to a certain extent. However, due to the complexity of acoustic wave and the limitation of Fourier transform [32], S and pseudo-Rayleigh waves are difficult to be completely suppressed. Figure 11c represents the relationship between ST wave energy and cave radius at different source distances. The energy of the ST wave gradually decreases with the increase of cave radius at different source distances.

3.6. The Effect of the Cave on the Scattered Wave

Figure 12a refers to the acoustic wave and the scattered wave. The wave propagation in the borehole of homogeneous formation is regarded as a direct wave, and the wave trains of cave formation in Figure 12a subtract the direct waves to separate the scattered waves. A Fourier transform was used to get the frequency spectrum of the scattered wave in Figure 12b. The peaks with larger amplitude are mainly S–S wave, and the dominant frequency is about 9–12 kHz. The peaks with smaller amplitude are mainly converted P-waves. We added the square of the value of each frequency to get the energy of the scattered wave. Taking the distance and radius of the cave as the abscissa and scattered wave energy as the ordinate, the relationship between the scattered wave energy and the cave was obtained in Figure 13.

The scattered wave recorded at 0.4 m below the cave (the source distance was 1.6 m) was selected. The cave with a diameter less than 0.25 times the wavelength of the S-wave was approximated as a point scatterer; otherwise the cave itself could form multiple-scattering [30]. In our model, a cave with a radius less than 2 cm could be regarded as a point-scatterer and could form multiple-scattering with a radius greater than 3 cm.

Figure 13a shows the relationship between the energy of scattered wave and the position of cave. When the cave can be regarded as a point-scatterer (R$\le$ 2 cm), the energy of the scattered wave has a negative correlation with the distance of the cave. The energy of scattered wave is very small when the distance is more than 24 cm and reaches the maximum when L = 2 cm if the cave itself can form multiple-scattering (R > 3 cm). The scattered wave energy gradually decreases when L > 2 cm (R > 3 cm).

Figure 13b shows the relationship between the energy of scattered wave and the radius of the cave. The energy of scattered wave raises with the increase of cave radius if R$\le$ 2 cm. If R > 3 cm, the energy of the scattered wave changes sharply when the cave is close to the borehole wall (L = 2 cm) and is relatively stable when L > 2 cm with the increase of cave radius.

4. Discussion

The numeric simulation of borehole acoustic waves in cave formations is a theoretical calculation; its ultimate purpose is to provide interpretation and reference for actual data. The results in our study show that borehole full-waves have response characteristics to the caves, which are mainly reflected in the energy and the first arrival of waves.

The porosity-frequency spectrum is calculated from the conductivity curve of electrical-imaging logging [33] and follows the following rules.

(1)

When the primary pores are well developed and the secondary pores are poorly developed, the porosity-distribution spectrum is a narrow single peak. The peaks move to the direction of small porosity;

(2)

When the primary pores are poorly developed and the secondary pores are well developed, the porosity-distribution spectrum is also a single peak. However, the peaks move to the direction of large porosity;

(3)

When the primary pores and secondary pores are well developed, but heterogeneous, the porosity-distribution spectrum shows multiple peaks. The peaks move to the direction of large porosity.

Cumulative porosity curves PS3, PS5, PS7…PS30 represent pore components with porosity values below 0.03, 0.05, 0.07…0.3 [34]. The porosity of electrical imaging can be obtained by filling the profile with different colors; the contribution of the different sizes of pores to the total porosity can be intuitively analyzed. According to the method in Figure 11b mentioned above, the energy of the ST wave is calculated. The acoustic full-wave was decomposed into several intrinsic mode function (IMF) components by empirical mode decomposition (EMD). P and S-wave were distributed in the highest frequency IMF1 component [35]. The energy of IMF1 component can be approximately regarded as the energy of the S-wave because the amplitude of the P-wave is very small. The slowness of the P-wave is the reciprocal of velocity. The first arrival delay of the P-wave indicates that the velocity is low, then the slowness is large.

Figure 14 shows the logging results of H-well in the igneous formation of Liaohe. Area a is basalt. Cores 1–2 show no alteration. The color of electrical-imaging logging is bright, indicating that the conductivity is low. The porosity-distribution spectrum shows a single peak moving forward, and the porosity of electrical imaging shows that small pores are mainly developed in the formation. The energy of the S and ST waves is high and stable. The density is large, and the porosity and slowness of the P-wave are small. All the above logging data indicate that this area is a tight formation. Cores 3–4 show altered basalt in area b. Dark spots appear in the electrical-imaging logging static image, which indicates that as conductivity increases, and caves develop. The porosity-distribution spectrum widens and moves in the direction of large porosity. It was found that there are mainly medium and large pores from the electrical imaging porosity. The energy of the S and ST waves was attenuated, and the slowness of the P-wave increased. This is in good agreement with the theoretical results of the cave formation simulated.

Figure 15 shows the effect of caves with different development degree on borehole acoustic waves. The attenuation of the S and ST waves becomes stronger, and the slowness of the P-wave becomes larger with an increasing number and size of caves.

The theoretical results of this study can be applied to practical data, but there are also some shortcomings.

(1)

The model simulated in our study is an ideal simplified model, which has a certain effect in describing the response characteristics of acoustic wave to a single cave. However, the number, shape and filling degree of caves in the actual formation are very complex. The author’s future work is to establish a 3D complex geological model and study the response characteristics of acoustic-logging more detailed and systematic;

(2)

The actual well-logging data can reflect the development degree of cave, but cannot accurately reflect the details of cave. It is a long way to go that how to use monopole acoustic-logging to detect reservoir caves precisely. More complex numeric simulation and comparison with actual data are needed;

(3)

A single GPU parallel-computing program is designed in this study. The memory of a single GPU may not meet the computing needs when the amount of data are too large, so we need to add support for multiple GPUs in the program. It is of great significance to apply CUDA parallel-computing of multiple GPUs to borehole acoustic numeric simulation.

5. Conclusions

The following important conclusions can be drawn:

• The investigation depth of monopole source was about 24 cm in our model;
• The cave cannot affect ST wave under the condition of low frequency. The position and size of cave affected the borehole acoustic wave at high frequency;
• Due to the superposition of scattered waves, the amplitude of the P-wave may increase or decrease. Theoretically, the first arrival of the P-wave can reflect the development degree of cave;
• The energy of the S-wave is always attenuated due to the propagation mechanism. S-wave is more sensitive to the presence of cave. The energy attenuation of the S-wave can reflect the size of the cave;
• The attenuation coefficient of the S-wave raises first, and then decreases due to energy conversion with the increase of source distance;
• For ST waves the amplitude decreases with the increase of cave radius;
• The excitation frequency may affect the relationship between the distance of the cave and P(S) waves;
• The energy of the scattered wave generally decreases with an increase of distance to the cave. It is worth noting that the scattered wave formed by large-scale near-borehole cave has no rule with radius and distance of the cave.