ABSTRACT

We developed a theoretical method of estimating the vertical and horizontal permeability of a transversely isotropic two-phase formation with a vertical symmetry principle axis (VTI) by using the attenuation characteristic of the dipole-flexural waves in a fluid-filled borehole. The attenuation value of flexural waves was extracted from synthetic data. The inversion was operated by the least-squares method. Two main conclusions were obtained based on the analysis of the results of inversion implemented on multiple sets of synthetic waveform data. First, the simultaneous estimation of vertical and horizontal permeability was feasible for the slow VTI formation, which means its shear-wave velocity is slower than the acoustic velocity of fluid in a borehole that is different from the previous understanding in some studies. Second, the fast formation corresponded to the opposite case from the slow formation, the precision of the estimated horizontal permeability relates to the accuracy of the vertical permeability, but the vertical permeability is difficult to estimate by the attenuation of the dipole-flexural waves. We added random noise with different intensities to synthetic data in order to simulate the errors that might exist in actual logging data. We adopted Prony’s method to extract attenuation to eliminate the effect of noise to a certain degree.

INTRODUCTION

Permeability is one of the most important parameters in geology to describe the production of an oil and gas reservoir. Accurately estimating permeability is one of the most basic tasks in oil and gas exploration and development, as well as one of the most important subjects of logging interpretation. Permeability has a relatively wide range of possible values and is difficult to measure. There is yet no single way with enough credibility to replace all other methods. The most prevalent current approach for determining permeability is via comprehensive interpretation and evaluation by several methods. So continuing to seek an effective new method to determine permeability is still a practically significant subject and a worthy direction for researchers to pursue.

It is generally believed that, compared to other measurement methods, the permeability value extracted from logging corresponds to the reservoir description scale and reflects in situ permeability better. Permeability is a mechanical parameter to describe fluid flow in two-phase media. The permeability measurement resulting from acoustic logging based on the principle of elastodynamics can be considered as a direct method. Based on the theoretical and field data analysis, the inference and prediction about the phase velocity and attenuation of the Stoneley waves in a borehole have correlation with formation permeability has been proposed by several researchers (Rosenbaum, 1974; White, 1983; Williams et al., 1984). Through a series of convincing laboratory experiments, and comparison with theory, Winkler et al. (1989) verify the validation of Biot’s model used to describe the dispersion and attenuation of the Stoneley waves in a borehole surrounded by permeable formation. The studies on estimating permeability quantitatively by applying Stoneley waves of acoustic logging were performed successively. Tang et al. (1991) realize rapid permeability estimation by the Stoneley wave in a simplified Biot-Rosenbaum model (Biot, 1956a; Rosenbaum, 1974) and took into consideration the influence of the logging tool size (Tang and Cheng, 1996). Wu et al. (1995) study the permeability inversion by the attenuation of the Stoneley wave based on the complete Biot-Rosenbaum model. Furthermore, the correct program for the impermeable intrinsic attenuation was proposed using the analysis of argillaceous by Wu et al. (1996). Liu and Johnson (1997) and Ticchelaar (1999) consider the effect of mud-cake on the Stoneley wave and achieved the quantitative inversion of permeability. Brie et al. (1998) realize the quantitative evaluation of formation permeability from Stoneley waves, taking into account the modification of the effect of mud and mud-cake properties. This methodology of estimating permeability by Stoneley waves was also used to process the actual field data (Pampuri et al., 1998).

In the last twenty years, the technology of acoustic logging has been developing rapidly that greatly expands the function of acoustic detection and provides a vast amount of measured data. Wang (1999) discovers that there is better correlation between the attenuation of the flexural wave and the permeability of two-phase medium formation, which is similar to the situation for Stoneley waves, and proposes the program of permeability inversion using the amplitude ratio of the flexural wave based on the study of forward modeling. The feasibility of their method has been verified by processing the synthetic signal and using practical examples, which offer a different way to get the permeability of formation from the Stoneley wave permeability. Modern devices used in acoustic logging are all combinations of monopole and cross-dipole array and other functional tools. So applying two sets of field data gotten from different acoustic logging methods to estimate identical reservoir parameters, such as shear velocity and permeability, should always be useful because they can confirm and complement each other.

However, the above results are all acquired in the case of isotropic formations. In fact, reservoirs and permeability are usually anisotropic. Therefore, it would be necessary in further study to estimate anisotropic permeability by acoustic logging. A transversely isotropic (TI) formation is one of the most common anisotropic formations because of the formation deposition and the directional distribution of the cracking. There are five independent elastic constants for a TI elastic solid, and there are three more constants that can be calculated for TI two-phase media. The TI media with the vertical symmetry axis (VTI) was studied earlier than other kinds of TI formations. Schmitt (1989) studies the problem of multipole acoustic fields in a fluid-filled borehole surrounded by a VTI two-phase formation.

Because of the demands of acoustic logging in a horizontal or a high-angle borehole surrounded by a TI formation, research on the numerical simulation of a 3D dipole transient acoustic field in TI media with arbitrary symmetry axis has been conducted (e.g., Liu and Sinha, 2000; He and Hu, 2009; Yan, 2011). The study on the forward model of acoustic logging in anisotropic media has been promoted, but there were few works on inversion parameters of anisotropic media. For example, there is still no way to quantitatively invert all the elastic modulus in either TI medium with the horizontal symmetry axis (HTI) or VTI medium. In addition, the detection and inversion of anisotropic permeability have barely been developed. One possible reason is that the dipole-flexural wave is not sensitive to vertical permeability (Schmitt, 1989). Recently, He et al. (2010) come to a similar conclusion through investigating the sensitivity coefficients of the flexural wave attenuation, and estimated the single horizontal permeability by simulated data based on a given vertical permeability. We also found that there is a big disparity of the sensitivity degree of vertical permeability to the flexural wave attenuation between the fast and slow formations (the slow formation means the velocity of a shear wave in the skeleton is slower than the acoustic velocity in the fluid of borehole VS<Vf the fast formation means the velocity of a shear wave in the skeleton is faster than the acoustic velocity in the fluid of the borehole VS>Vf). The sensitivity of flexural wave attenuation in a slow VTI two-phase formation to vertical permeability is significantly greater than it is in fast VTI two-phase formation. Furthermore, the slow formation is actually a kind of important common reservoir also. Therefore, the synchronous inversion of vertical and horizontal permeability for VTI two-phase formation is still possible and should be explored continuously, at least for the slow formation.

Detecting formation anisotropy is an important task of cross-dipole acoustic logging. Developing the functions of the cross-dipole acoustic logging from detecting only velocity anisotropy to determining permeability anisotropy is a natural and logical advance in methodology. In continuously cross-dipole acoustic logging, many kinds of anisotropic formations, including HTI, VTI, or even more complicated anisotropic formations, can all be encountered. Using cross-dipole acoustic logging to realize the inversion of anisotropic permeability is a more reasonable choice in practical application also. This does not exclude the effect of the monopole acoustic logging. In fact, the completion of the inversion of anisotropic formation parameters usually needs the combination of two methods.

In accordance with the discussion concerning the possibility of the simultaneous inversion of vertical and horizontal permeability as outlined above, we focus our study on the inversion of anisotropic permeability of VTI two-phase media by dipole-flexural waves, which is also a beginning for solving more complicated problems of anisotropic permeability.

As an important stage before practical application, this paper only describes theoretical methods and attempts to invert vertical and horizontal permeability simultaneously and emphasizes for slow formation. Several works realize permeability inversion by Stoneley waves in a borehole (Tang et al., 1991; Wu et al., 1995, 1996; Tang and Cheng, 1996; Brie, 1998), which were all based on the traditional model studied to a higher frequency range (Biot, 1956b; Johnson, 1987), and improve it being applicable to anisotropic (VTI) two-phase media (Schmitt, 1989; Zhang 1995). To simulate the influence of possible errors in practical logging data, we added random noise of different intensities to them. Finally, the vertical and horizontal permeability inversion for fast formations using synthetic data is also conducted and discussed.

METHOD

Borehole configuration

The borehole is assumed to be a fluid-filled cylinder with the radius a. The density and acoustic speed of the fluid in the borehole are ρf and Vf, respectively. The medium outside the borehole is the transversely isotropic two-phase medium, whose symmetric axis is parallel to the borehole axis. The dipole source and receivers, which are used to excite and receive the acoustic field, are all located on the borehole axis (Figure 1).

The displacement potentials

We adopt the cylindrical coordinate system (r, θ, and z) centered on the dipole source and oriented along the borehole axis. Because the excitation frequency of the acoustic logging tool is obviously higher than seismic exploration frequency and covers a wider frequency range, Biot’s (1956b) high-frequency theory is adopted (Johnson, 1987). The dynamic differential equations in frequency domain are as follows:  
{·σ=2t2(ρ11·U+ρ12·V)+b·FJ(ω)·t(UV),σ=2t2(ρ12·U+ρ22·V)b·FJ(ω)·t(UV),
(1)
where σ is the stress tensor for the skeleton and σ is the effective fluid stress for saturated fluid (Schmitt, 1989; Zhang, 1995). Here, ρ11, ρ12, and ρ22 are the mass coefficients tensor, and these coefficients are considered when the borehole is not uniform:  
{ρ11=(1ϕ)ρsI(IE)ϕρfρ12=(IE)ϕρfρ22=Eϕρf,
(2)
where I is the unit tensor and ρs and ρf represent the skeleton density and fluid density. Where ϕ is the porosity of media and E is the intrinsic tortuosity tensor, which is defined by 
E=(EH000EH000EV).
(3)
Symbol b is the second-order tensor, which is defined by 
b=(ηϕk),
(4)
where η is the intrinsic viscosity and k is the intrinsic anisotropic permeability tensor defined by 
k=(kH000kH000kV),
(5)
where FJ(ω) is the viscous corrected vector factor, which should be considered when the viscous skin depth becomes comparable to the relevant pore size (Biot, 1956b; Johnson, 1987) that is taken for the expression by Johnson (1987): 
FJ(ω)=(FJH(ω)000FJH(ω)000FJV(ω)){FJH(ω)=14iEHkHρfωmηϕFJV(ω)=14iEVkVρfωmηϕ,
(6)
where m is a nondimensional parameter, and usually takes a value from 8 to 12 (Johnson, 1987).
With the viscous correction factor, the viscosity is frequency dependent and so is permeability. Here, U and V represent displacements of the skeleton and fluid, respectively. Then, we define γ as follows: 
γmmH=ρmmH+iηϕ2ωkHFJH(ω),γmmV=ρmmV+iηϕ2ωkVFJV(ω)(m=1,2),γ12H=γ21H=ρ12Hiηϕ2ωkHFJH(ω),γ12V=γ21V=ρ12Viηϕ2ωkVFJH(ω),
(7)
where the superscripts H and V correspond to the horizontal and vertical components of the parameters.
With the time dependence factor eiωt, equation 1 can be expressed in the frequency domain as 
·σ=ω2(γ11HUr+γ12HVrγ11HUθ+γ12HVθγ11VUz+γ12VVz)σ=ω2(γ12HUr+γ22HVrγ12HUθ+γ22HVθγ12VUz+γ22VVz).
(8)
The relationships of stress and strain in the two-phase VTI formation are as follows: 
(σrrσθθσzzσθzσrzσrθσ)=(AA2NF000MA2NAF000MFFC000Q000L0000000L0000000N0MMQ000R)(erreθθezz2eθz2erz2erθε),
(9)
where A, N, C, F, and L are the elastic modulus of the TI formation, M and Q are for solid-fluid coupling, and R is for the fluid-phase. Here, eij and ε represent the strain components of the solid and fluid, respectively. The M, Q, and R are calculated by  
{M=ϕ(α1ϕ)[(αϕ)/Ks+ϕ/Kf]Q=ϕ(α3ϕ)[(αϕ)/Ks+ϕ/Kf]R=ϕ2[(αϕ)/Ks+ϕ/Kf],where{α=(2α1+α3)3α1=1(2A2N+F)3Ksα3=1(C+2F)3Ks,
(10)
where Ks and Kf are bulk compressional modulus of grains and fluid and Kf=Vf2ρf, respectively.
The displacement potentials are defined by the formulas below: 
U=φ+×(χez)+××(ψez),V=φ+×(χez)+××(ψez),
(11)
where φ and φ are compressional potentials, χ and χ are horizontal polarized shear potentials (SH-wave), and ψ and ψ are vertical polarized shear potentials (SV-wave).

Wave potentials

Based on the analysis above, we get the solutions of displacement potential equations by a set of derivations, (given in Appendix A), as  
χ=(νr0)n2nn!χ0Kn(αSHr)sinn(θθ0),
(12)
 
χ=(νr0)n2nn!a0χ0Kn(αSHr)sinn(θθ0),
(13)
 
φ=(νr0)n2nn![a1g1Kn(α1r)+g2Kn(α2r)+g3Kn(α3r)]cosn(θθ0),
(14)
 
φ=(νr0)n2nn![b1g1Kn(α1r)+b2g2Kn(α2r)+b3g3Kn(α3r)]cosn(θθ0),
(15)
 
ψ=(νr0)n2nn![g1Kn(α1r)+a2g2Kn(α2r)+a3g3Kn(α3r)]cosn(θθ0),
(16)
 
ψ=(νr0)n2nn![c1g1Kn(α1r)+c2g2Kn(α2r)+c3g3Kn(α3r)]cosn(θθ0).
(17)

Equations 12–13 are the solutions of SH-wave, and equations from 14 to 17 are the solutions of P–SV-wave (P- and SV-wave are coupled strongly in the transversely isotropic media). The In and Kn are the nth modified Bessel functions of the first and second kind. We use the symbol ν to represent the radial wavenumber in fluid (ν=kz2kf2), kz represents the axial wavenumber, and kf=ω/Vf. For the SH-wave (equations 12 and 13), χ0 is the weighting coefficient, αSH is the radial wavenumber, and a0=γ12H/γ22H (the meaning of γ is given in Appendix A). For the P–SV-wave (equations 14–17), gi and αi(i=1,2,3) express weighting coefficients and radial wavenumber of quasi-P1(i=1), quasi-P2(i=2), and quasi-SV wave (i=3), respectively. The ai, bi, and ci(i=1,2,3) are unknown coefficients, which can be solved by substituting equations from 14 to 17 to the boundary conditions. The order of the multipole source is denoted by n, the flexural wave is excited by the dipole source, which corresponds to the condition of n=1, which is discussed below, so the subscript n is omitted. Based on equations 12–17, the field equations A-12, A-23–A-28 can be derived, and the derivation is given in Appendix A.

Then we can get the representations of displacement and stress, which are shown in Appendix A.

The representation of displacement potential of the fluid inside the borehole is given by 
ϕf(r,θ,kz,ω)=νr02[2K1(νr)+AI1(νr)]cos(θθ0).
(18)
Here, r0 and θ0 are the separation and azimuth angle of the dipole source. Based on equation 18, the associated radial displacement and the stress of fluid-phase can be solved (equations A-29 and A-30), the specific process is given in Appendix A.

Wave propagation

In this paper, we consider the interface between the fluid-filled borehole and the VTI two-phase formation as a permeable wall only, so the boundary conditions at the position r=a of the cylindrical coordinate system are 
UrI=UrII+ϕ(VrIIUrII),pfI=TrrII,pfI=σIIϕ,0=TrzII,0=TrθII,
(19)
where the superscript I and II represent the media inside and outside borehole, respectively. Here, T represents the total stress, which is defined by Tij=σij+σδij. We substitute the acoustic-field equations (from A-10–A-12 and A-23–A-30) into equation 19, and we obtain the following equation: 
(m11m12m13m14m15m21m22m23m24m25m31m32m33m34m350m42m43m44m450m52m53m54m55)(A1g2g3g1χ0)=(b1b2b300),
(20)
where the elements mij and b1,b2,b3 of the matrix above are given in Appendix B. The reflection coefficient A1 can be represented as  
A1(kz,ω)=N(kz,ω)D(kz,w);
(21)
here, 
D(kz,ω)=(b1m12m13m14m15b2m22m23m24m25b3m32m33m34m350m42m43m44m450m52m53m54m55),N(kz,ω)=(m11m12m13m14m15m21m22m23m24m25m31m32m33m34m350m42m43m44m450m52m53m54m55).
(22)
Then, the acoustic field excited by a dipole source in the borehole in frequency–wavenumber domain can be written as 
P(kz,ω)=kfr0F(ω)2[A1(kz,ω)I1(kfr)+4K1(kfr)]cos(θθ0);
(23)
here, F(ω) is the frequency spectrum of the source excitation pulse.

Inversion method

The flexural wave is a type of guided wave whose information is contained in the reflection coefficient. The poles of equation 21 correspond to all the modes of flexural waves. The amplitude of the flexural wave in the spatial-frequency domain can be calculated by the poles’ residue and can be written as 
Pfl(z,ω)=2πiF(ω)(kfr02)I(kfr)cos(θθ0)Res(N(kz,ω)D(kz,ω))eikzz|kz=kzn.
(24)

The flexural wave has attenuation characteristics because the VTI two-phase formation outside the borehole is dissipative (White, 1983; Tang and Cheng, 2004). Thus, the wavenumber of flexural wave is complex, which is expressed by kzfl.

Let 
kzfl=k+Riα(ω),
(25)
where kR and α represent the propagation and attenuation of the flexural wave, respectively. It is conceivable that the attenuation can be extracted by amplitude ratio of the flexural wave in theoretical study. However, the practical logging data are more complicated; this method seems too be simple to extract accurate attenuation. To approach a practical application, we use Prony’s method to extract attenuation of the flexural waves. Prony’s method has become an effective way to estimate power spectral density (Kay and Marple, 1981). Moreover, many researchers applied Prony’s method to estimate slowness dispersion and attenuation from arrays of sonic logging waveforms (Parks et al., 1982, 1983; Lang, 1987; Ma et al., 2010).

First, we analyze the effect of the intrinsic anisotropic permeability of VTI two-phase media on the flexural wave attenuation by numerical simulation.

The elastic modulus of the VTI two-phase media are given in Table 1, these parameters are quoted from Schmitt (1989). We take three and four for the intrinsic horizontal and vertical tortuosity, and nine for m. Theoretically, the value of tortuosity relates to the value of permeability, but that effect is so little that can be ignored. Therefore, we always take these values in the calculation of different permeabilities. For reflecting the features of slow and fast formation and their anisotropy directly, the parts of characteristic wave velocities in the VTI two-phase medium are shown in Table 2; they include fast and slow P-waves, SV-, and SH-waves, respectively. These velocities are in the condition of 6-kHz frequency with the intrinsic vertical and horizontal permeability — both 300 mD. Both the fluid-phase in the two-phase formation and the fluid in the borehole are assumed to be water. The borehole radius and the dipole source separation were 0.1 and 0.01 m, respectively. In VTI two-phase formation, the permeability is related to the directions, so we make kH and kV represent the intrinsic horizontal and vertical permeability, respectively.

Figure 2 displays the relationship between the attenuation α(ω) of the flexural wave and different permeability values in the frequency domain. The intrinsic permeability values kH and kV are changed from 10 to 1000 mD in the slow and fast VTI two-phase formations. It shows from Figure 2 that the attenuation α(ω) increases along with the frequency. The α(ω) relates to intrinsic permeability values kH and kV. The effect of the permeability on α(ω) in the high-frequency area is greater than that in the low-frequency area. It could be found that the effect of kH on α(ω) is significantly greater than that of kV on α(ω), which is consistent to the conclusion that the attenuation is not sensitive to kV. Figure 2 also shows that there is a significant difference in the sensitivity between the slow and fast formations. The attenuation in the slow formation is much more sensitive to the vertical permeability than that in the fast formation. This is a favorable conclusion for estimating the vertical and horizontal permeability simultaneously in the slow formation. In this paper, we focus our attention on the intrinsic anisotropic permeability inversion in the slow formation that covers some of the known oil reservoirs. Meanwhile, we also study the intrinsic anisotropic permeability inversion by the dipole-flexural wave in fast VTI two-phase media as a comparative work.

The inversion analysis is then carried out. The attenuation of the flexural wave, which is extracted from the time-domain full waveforms by Prony’s method, is named measured value and denoted as D(ω,kH,kV), where kH and kV are the real intrinsic horizontal and vertical permeability values in the formation. The attenuation α(ω) calculated by equation 25 expressed by G(ω,kH,kV), in which kH and kV are the estimated intrinsic permeability values. Other parameters of the formation are given values. Then, D and G are substituted in the variance function: 
Q=i=1N[G(ωi,kH,kV)D(ωi,kH,kV)]2,
(26)
as the target function, where N is the number of the frequency points in selection. The least-squares method of double unknown quantities is used to invert kH and kV.

Feasibility investigation

To ensure the smoothness of inversion and the adaptability of the method, the variation function Q was studied for different values of estimated intrinsic permeability (kH and kV). The vertical and horizontal permeability values of slow VTI two-phase formation are investigated, and they are assumed to be 400 and 500 mD, respectively. We selected 80 points in the frequency domain from 4 to 8 kHz (N=80), the range of kH is from 200 to 600 mD and that of kV is from 200 to 800 mD. The variation function Q is calculated when kH and kV are changed, and the result is shown in Figure 3.

From Figure 3a, we can see that the value of Q has many poles when kH and kV changes, and reaches the minimum value when kH=kH and kV=kV. It is clearer in its plane projection image in Figure 3b, the values of kH and kV in the minimum pole correspond to 400 and 500 mD, respectively. This is equal to the true value in the formation, which means the variation function Q reaches the global minimum when the estimated permeability is close to the true value. Through changing initial values of kH and kV to test the convergence of the inversion method, for different true values kH and kV, trying many times, it always converges smoothly to the true values, even though there are some local minimum in the neighborhood of the global minimum as shown in Figure 3. A possible reason is that the variation function Q of local minimums is too large to satisfy the least-squares method. In conclusion, the convergence of the inversion method is smooth and not sensitive to the initial values of kH and kV in quite a wide range. It is feasible to estimate permeability for the theoretical model by the least-squares method.

Therefore, the least-squares method is suitably used to invert the vertical and horizontal permeability simultaneously in slow and fast VTI two-phase media. The measured value of the attenuation can be extracted from the flexural wave of dipole full waveforms by Prony’s method. We observe that the flexural wave calculated by pole residue is basically the same as that extracted from the full waveforms, so the flexural wave calculated by equation 24 can be used as actual waveform as a theory study. The attenuation extracted from the full waveforms above could be seen as measured data in inversion study. Our aim is to numerically investigate the feasibility, adaptability, and accuracy of estimating the anisotropic permeability by the flexural wave attenuation in VTI two-phase formation.

RESULTS OF NUMERICAL INVESTIGATION

In this section, the simultaneous inversion of vertical and horizontal permeabilities are performed on a series of synthetic waveforms data, and the effect of noise on the inversion results is also discussed numerically.

First, we estimate permeability by attenuations extracted from synthetic data of dipole-flexural waves. The inversion is implemented by the least-squares method, and its convergence is very smooth for any initial value.

Table 3 shows the inversion results of different vertical and horizontal permeability of the slow VTI two-phase formation whose elastic modulus are given in Table 1. Here, kH and kV are the true intrinsic permeability values of the simulation reservoir, and kH and kV are the inversion results. Meanwhile the relative errors of them are also shown in the third and sixth rows of Table 3. We infer that horizontal and vertical permeabilities of slow porous formations can be estimated simultaneously based on our previous investigation. The curves of different horizontal permeabilities clearly separate from each other (Figure 2a) and that of the different vertical permeabilities separate a little (Figure 2b). It can be seen that the results coincident with our deduction, the relative errors of inversed horizontal permeability are all less than 0.3%, whereas that of inversed vertical permeability are under 2% too. The inversing precision is very high, and the convergence of inversion is quite good. We notice that when the horizontal permeability is smaller than the vertical permeability, the relative error of the inversion horizontal permeability is usually large, such as the sixth row and eighth row of Table 3. Therefore, we preliminarily infer that the error of the inversion result may increase when the vertical permeability is much larger than the horizontal permeability. That means the inversion results of vertical and horizontal permeability are influencing each other; if one of them is inaccurate, the other would be far away from the true value, too.

The inversion discussed above is idealized and indicates that it is possible to invert vertical and horizontal permeability simultaneously in theory. However, at this point, we cannot affirm yet that this method is feasible in practical application. The actual logging data are much more complicated than the synthetic data because the theoretical model is always more simple than the actual formation, and there is a higher probability of systematic and accidental errors. To reflect the error of the actual data, we investigate the inversion results on the synthetic data with random noise.

We add 5% random noise to waveforms of the five receivers with equal spacing (0.5 m) in the same formation. We adopt Prony’s method to extract attenuation of flexural waves to eliminate the effect of noise. Generally, it is believed that, compared to other acoustical characteristics, velocity, for example, the accuracy of extracted attenuation is much more affected by noise. The extracted attenuation is unstable in low- and high-frequency areas, except around the central frequency range. We abandon these unstable data and only select the middle frequency area (5–6 kHz). It definitely promotes the inversion precision. We then invert several groups of vertical and horizontal permeability in the slow two-phase formation by the method discussed above. The results are shown in Table 4. We can see that when the true horizontal and vertical permeability is relatively small (50 mD), the relative error of the inversion horizontal permeability result is a little large, but the absolute error of which is acceptable. On the contrary, the inversion result of vertical permeability on the low permeability condition is far from true value. Low permeability causes small attenuation, which means the attenuation barely changes, whereas vertical permeability changes. As the true permeability increases, the relative errors of inversion results decrease rapidly. For horizontal permeability, the relative errors are all smaller than 10%, most of them are less than 4%. The precision can be accepted in practical applications. For vertical permeability, except the low permeable situation, the relative errors are acceptable too. Most of them are less than 10%. We also notice that, corresponding to Table 3, the inversion results’ error of the sixth and eighth rows are relatively large as well. This is an interesting phenomenon, and it partly proves our conclusion that the precision of inversion results of vertical and horizontal permeability interact somehow. That illustrates the importance of estimating horizontal and vertical permeability simultaneously.

Furthermore, the receiver array with eight or 12 receivers is usually applied in a practical logging system, and the spacing is 0.15 m, as usual. We adopt a receiver array with eight receivers to invert the permeability, to approach the actual application. With spacing of 0.15 m instead of 0.5 m previously, the attenuation will be smaller and more difficult to extract. The inversion results are shown in Table 5. Compared with the results of five receivers, the accuracy of inversions are not as good as Table 3, especially for vertical permeability. In addition, as in Table 3, the relative error of inversion horizontal permeability is large, whereas that of inversion vertical permeability is large. We infer that the shorter spacing of receivers causes the larger inversion error of the results.

Finally, the vertical and horizontal permeability of a fast VTI two-phase formation are also inversed by theoretical simulation data, and the results are shown in Table 6. In Figure 2c, we can see that the attenuation curves of different horizontal permeabilities separate clearly, but that curves of different vertical permeabilities almost superpose in Figure 2d. The results are in accordance with the expectation that the inversion precision of horizontal permeability is high. It is notable that most of the vertical permeability results are credible also, these results greatly exceed our expectation. Though the method still needs improvement before practical application because the practical waveforms are more complicated than the synthetic waveforms, we could still get enlightening results. From Table 6, we can see that the errors of some horizontal permeability inversion are larger than others’, whereas errors of the corresponding vertical permeability are large too, and the relative error is high when vertical permeability value is larger than horizontal permeability (kV>kH). This illustrates that the precision of the inversion of horizontal permeability strongly relates to the vertical permeability. This conclusion is suitable for slow VTI two-phase formation too. In practical inversion of VTI two-phase formations, a priori value of vertical permeability would affect directly on the precision of inversion results of horizontal permeability when the vertical and horizontal permeability could not invert simultaneously.

DISCUSSION

The estimation of anisotropic permeability of a reservoir is an interesting and challenging subject. Because the inversion of vertical permeability of the VTI two-phase formation using the dispersion and attenuation of dipole-flexural waves is considered impossible according to previous understandings (Schmitt, 1989; He and Hu 2009; He et al., 2010), the study in this paper is still focused on the simultaneous estimation of vertical and horizontal permeability of VTI two-phase formations in theory.

In the next step of the actual data processing, we have to solve the influence of a variety of factors on the permeability inversion. We think that the influence factors mainly include two aspects. One is the value assignment of other parameters except unknown permeabilities; their value can be obtained by petrophysics measurements and other wireline logs as well as both monopole and dipole acoustic logging, and the analysis of the sensitive coefficient to some quantity can usually simplify treatment process. The other one is the change or damage of the physical environment, among which the influence of the mud cake and impairment layer of a borehole is a particularly important and perplexing problem. The flexural waves with horizontal polarization propagating along the skeleton of a poroelastic formation are different from the character of Stoneley waves; therefore, a new study on this topic should be developed. Actually, these problems can only be solved combining with the actual data processing, so making a feasible strategy (technical route) of the inversion would be necessary.

CONCLUSION

In this paper, we carry out a theoretical and numerical investigation of the simultaneous estimation of the intrinsic vertical and horizontal permeability for VTI two-phase (Biot) formations using the attenuation characteristic of the dipole-flexural waves in a fluid-filled borehole. The attenuation value of flexural waves is extracted by Prony’s method from synthetic waveforms data of dipole-flexural waves in frequency domain. Based on the analysis of the extreme value distribution of variance functions structured by model function and measuring data, the inversion of vertical and horizontal permeability is implemented using the least-squares method; the convergence of the inversion for multiple sets of parameters are all smooth.

For the slow VTI two-phase formation, the precision of inversion results for multiple sets of permeability are all higher, which means it is feasible to estimate vertical and horizontal permeability simultaneously for a slow VTI two-phase formation, at least in theory. The influence of error that is possibly presented in field data on the inversion precision is investigated numerically by adding random noise to synthetic data. Noise largely effects the accuracy of extracted attenuation, especially on the two poles of the frequency band. Thus, we abandon the data of the low- and high-frequency areas to eliminate the noise effect.

For the fast VTI two-phase formation, the inversion of vertical and horizontal permeability is tentatively carried out also for synthetic data. It is found that the inversion precision of horizontal permeability strongly relates to the precision of vertical permeability. This is a new understanding, which means that although the attenuation of dipole-flexural waves is not very sensitive to vertical permeability of a fast formation, it is not yet expected that a precise result of horizontal permeability can be obtained by single parameter inversion, no matter how much the true value of vertical permeability is.

ACKNOWLEDGMENTS

This work is supported by the National Nature Science Foundation of China (grant nos. 11134011 and 11174322). We thank G. Harrison for his earnest help in corrections for language usage.

THE DERIVATION OF FIELD EQUATIONS

Six displacement potential equations can be obtained by equations 8–11: 
N2χ+(LN)2χz2+ω2(γ11Hχ+γ12Hχ)=0,
(A-1)
 
γ12Hχ+γ22Hχ=0,
(A-2)
 
A2φ+(F+2LA)2φz2+M2φ+ω2(γ11Hφ+γ12Hφ)+z[(ALF)2ψ+(F+2LA)ψ2z2+ω2(γ11Hψ+γ12Hψ)]=0,
(A-3)
 
z[(F+2L)2φ+(CF2L)2φz2+Q2φ+ω2(γ11Vφ+γ12Vφ)]+(2z22)[L2ψ+(CF2L)2ψz2+ω2(γ11Vψ+γ12Vψ)]=0,
(A-4)
 
M2φ+R2φ+(QM)2φz2++ω2(γ21Hφ+γ22Hφ)+z[(QM)(2z22)ψ+ω2(γ21Hψ+γ22Hψ)]=0,
(A-5)
 
z[(γ21Hγ21V)φ+(γ22Hγ22V)φ]+2z2(γ21Hψ+γ22Hψ)+(22z2)(γ21Vψ+γ22Vψ)=0.
(A-6)

The displacement potential equations of the SH-wave are equations A-1 and A-2, so we can easily get equations 12 and 13 from equations A-1 and A-2.

The P–SV-wave displacement potential function are equations A-3–A-6; we can take the trial solution of them in the form below:  
φ=R1Kn(αr)cosn(θθ0),φ=R2Kn(αr)cosn(θθ0),ψ=R3Kn(αr)cosn(θθ0),ψ=R4Kn(αr)cosn(θθ0).
(A-7)
Equations 11–14 can be solved as follows:  
[Aα2(F+2L)kz2+ω2γ11h]R1+[M(α2kz2)+ω2γ12h]R2+ikz[(AFL)α2Lkz2+ω2γ11h]R3+ikzω2γ12hR4=0,
(A-8)
 
ikz[(F+2L)α2Ckz2+ω2γ11v]R1+ikz[Qα2Qkz2+ω2γ12v]R2α2[Lα2(F+LC)kz2+ω2γ11v]R3α2ω2γ12vR4=0,
(A-9)
 
(Mα2Qkz2+ω2γ12h)R1+[R(α2kz2)+ω2γ22h]R2+ikz[(MQ)α2+ω2γ12h]R3+ikzω2γ22hR4=0,
(A-10)
 
ikz(γ12hγ12v)R1+ikz(γ22hγ22v)R2+(α2γ12vkz2γ12h)R3+(α2γ22vkz2γ22h)R4=0.
(A-11)
The coefficient determinant of equations A-8–A-11 must be zero to ensure R1, R2, R3, and R4 are nonzero. Then we get 
d1α6+(d2ω2+d3kz2)α4+(d4ω4+d5kz4+d6ω2kz2)α2+d7ω6+d8ω4kz2+d9ω2kz4+d10kz6=0,
(A-12)
where di(i=1,,10) are in the form below:  
d1=Lγ22V[AR(M)2],
(A-13)
 
d2=[ARM2][γ11Vγ22V(γ12V)2]+Lγ22V[Aγ22H+Rγ11H2Mγ12H],
(A-14)
 
d3=Lγ22H(M2AR)+γ22V[2L(FRMQ)+F(FR2MQ)+C(M2AR)+AQ2],
(A-15)
 
d4=[Aγ22H+Rγ11H2Mγ12H][γ11Vγ22V(γ12V)2]+Lγ22V[γ22Hγ11H(γ12H)2],
(A-16)
 
d5=Lγ22V(RCQ2)+γ22H[2L(MQFR)+F(2MQFR)+C(ARM2)AQ2],
(A-17)
 
d6=RL[γ11Hγ22H(γ12H)2]+γ11Hγ22V(Q2RC)+γ11Vγ22H(M2AR)+γ22Hγ22V(2FL+F2AC)+2γ12Hγ22V[MCQ(F+L)]+2γ12Vγ22H[QAM(F+L)]+2γ12Hγ12V[R(F+L)MQ],
(A-18)
 
d7=[γ11Hγ22H(γ12H)2][γ11Vγ22V(γ12V)2],
(A-19)
 
d8=(2Qγ12VCγ22VRγ11V)[γ11Hγ22H(γ12H)2]Lγ22H[γ11Vγ22V(γ12V)2],
(A-20)
 
d9=[RC(Q)2][γ11Hγ22H(γ12H)2]Lγ22H(Cγ22V+Rγ11V2Qγ12V),
(A-21)
 
d10=Lγ22H[(Q)2CR].
(A-22)
There are always three solutions of α in equation A-12 corresponding to α1,α2 and α3, respectively. For a given αi(i=1,2,3), only one independent value of R1,R2,R3,R4 exists. Meanwhile α1,α2 and α3 correspond to quasi-P1, quasi-P2, and quasi-SV, respectively. Because P1 and P2 waves are coupled with SV-waves strongly in TI media, we can gain the equations from 13 to 17.
Then, based on the six potential equations (equations 12–17), the displacement and stress of the formation outside the borehole are in the following form: 
Ur=νr02[(a1+ikz)α1K1(α1r)g1+i=23(1+ikzaj)αjK1(αjr)gj+nrχ0K1(αSHr)]cos(θθ0),
(A-23)
 
Vr=νr02[(b1+ikzc1)α1K1(α1r)g1+i=23(bj+ikzcj)αjK1(αjr)gj+1ra0χ0K1(αSHr)]cos(θθ0),
(A-24)
 
Trz=νr02L{(2ikza1α12kz2)α1K1(α1r)g1+i=23[2ikz(αj2+kz2)aj]αjK1(αjr)gj+irkzK1(αSHr)χ0}cos(θθ0),
(A-25)
 
Trr=νr02{[[(A+M)(a1+ikz)α12+ikz(F+Q)(ikza1α12)+(M+R)(α12kz2)b1+2r2N(1+ikzaj)n(n+1)]K1(αjr)+2rN(1+ikzaj)αjK2(αjr)]+2rN[n1rK1(αSHr)αSHK2(αSHr)]χ0}cos(θθ0),
(A-26)
 
Trθ=Nνr02{2r(a1+ikz)α1K2(α1r)d12rj=23(1+ikzαj)αjK2(αjr)dj+[kSH2Kn(αSHr)+2rαSHK2(αSHr)]χ0}sin(θθ0),
(A-27)
 
σ=vr02{[(α12kz2)(a1M3+M6b1)+ikz(M2M3)(ikza1α12)+ikz(M5M6)(ikzb1c1α12)]K1(α1r)g1+j=23[(αj2kz2)(M3+M6bj)+ikz(M2M3)(ikzαj2aj)+(M5M6)ikz(ikzbjαj2cj)]K1(αjr)gj}cos(θθ0).
(A-28)

Here, T represents the total stress tensor, which is defined by Tij=σij+σδij.

For the fluid field in the borehole, the associated radial displacement Urf and the stress Pf can be derived by the potential function (equation 18). Given that Pf=ω2ρfϕf and Urf=ϕf, we determine the representations of fluid displacement and stress inside the borehole as follows:  
Urf=ϕfr=νr02[2νK1(νr)+νAI1(νr)]cos(θθ0),
(A-29)
 
Pf=ω2ρfνr02[2K1(νr)+AI1(νr)]cos(θθ0).
(A-30)

THE ELEMENTS OF THE MATRICES IN EQUATION (20)

 
m11=νIn(νa),
(B-1)
 
m12=[(1φ)(1+ikza2)+φ(b2+ikzc2)]α2Kn(α2a),
(B-2)
 
m13=[(1φ)(1+ikza3)+φ(b3+ikzc3)]α3Kn(α3a),
(B-3)
 
m14=[(1φ)(a1+ikz)+φ(b1+ikzc1)]α1Kn(α1a),
(B-4)
 
m15=(1φ+φa0)naKn(αSHa),
(B-5)
 
m21=ρfω2In(νa),
(B-6)
 
m22=[(A+M)(1+ikzc2)α22+ikz(F+Q)(ikzα22c2)+(M+R)(α22kz2)a2+2a2N(1+ikzc2)n(n+1)]Kn(α2a)+2aN(1+ikzc2)α2Kn+1(α2a),
(B-7)
 
m23=[(A+M)(1+ikzc3)α32+ikz(F+Q)(ikzα32c3)+(M+R)(α32kz2)a3+2a2N(1+ikzc3)n(n+1)]Kn(α3a)+2aN(1+ikzc3)α3Kn+1(α3a),
(B-8)
 
m24=[(A+M)(a1+ikz)α12+ikz(F+Q)(ikza1α12)+(M+R)(α12kz2)b1+2a2N(a1+ikz)n(n1)]Kn(α1a)+2aN(a1+ikz)α1Kn+1(α1a),
(B-9)
 
m25=2aNa[n1aKn(αSHa)αSHKn+1(αSHa)],
(B-10)
 
m31=ρfω2In(νa),
(B-11)
 
m32=1φ[(α22kz2)(M+Ra2)+ikz(QM)(ikzα22c2)]Kn(α2a),
(B-12)
 
m33=1φ[(αj32kz2)(M+Ra3)+ikz(QM)(ikzα32c3)]Kn(α3a),
(B-13)
 
m34=1φ[(α12kz2)(a1M+Rb1)+ikz(QM)(ikza1α12)]Kn(α1a),
(B-14)
 
m42=L[2ikz(α22+kz2)c2]α2Kn(α2a),
(B-15)
 
m43=L[2ikz(α32+kz2)c3]α3Kn(α3a),
(B-16)
 
m44=L[2ikza1α12kz2]α1Kn(α1a),
(B-17)
 
m45=LiakznKn(αSHa),
(B-18)
 
m52=2Nan(1+ikzc3)[n1aKn(α2a)α2Kn+1(α2a)],
(B-19)
 
m53=2Nan(1+ikzc3)[n1aKn(α3a)α3Kn+1(α3a)],
(B-20)
 
m55=N[αSH2+2na2(n1)]Kn(αSHa)+2NaαSHKn+1(αSHa),
(B-21)
 
b1=εnnaKn(νa)εnνKn+1(νa),
(B-22)
 
b2=b3=εnρfω2Kn(νa).
(B-23)
Freely available online through the SEG open-access option.