## ABSTRACT

Near a borehole, stress concentration effects may cause a complex spatial variation of elastic anisotropy. Stress-induced sonic anisotropy results when moduli and velocities are stress dependent and the state of stress is nonhydrostatic. For such cases, the appropriate material model is that of an elastic orthorhombic medium with material axes aligned with the principal axes of stress. We simulate the dispersion characteristics of quadrupole waves propagating along a borehole in a stress-sensitive formation. To include sensitivity to stress in the analysis, we derive a finite-element modeling approach which includes the stress-velocity relationship, the perturbed stress and velocity field near the borehole, and solutions to the full coupled anisotropic and 3D wave equation. We characterized the anisotropy around the borehole for different in situ stress conditions, and we derived the effects on quadrupole dispersion characteristics. Of particular interest was the phase velocity of the quadrupole wave, which at low frequencies provides an estimate of formation shear-wave velocity. Results indicated that the quadrupole mode dispersion was relatively insensitive to the hoop stress effects around a borehole giving greater confidence to shear-wave estimates made from these modes. The proposed ring-element formulation, which includes Fourier expansions of the field variables in the direction of the circumference, is well suited to this analysis and can be adapted to include other effects, such as a finite-length borehole.

## INTRODUCTION

Quadrupole modes are guided waves that propagate along a borehole with the azimuthal acoustic profile $cos(2\theta )$ in the propagating medium. In isotropic and homogeneous formations, the velocity of quadrupole waves moves asymptotically to the Scholte wave speed at high frequencies and to the formation shear speed at low frequencies. There is practically no energy directed into this wave at frequencies below the cut-off frequency velocity, where the notation of a cut-off frequency is used to refer to the frequency at which the phase velocity of the mode is equal to the formation shear velocity. Due to the asymptotic behavior, quadrupole acoustic logging can be used for estimating in situ shear-wave velocities in real formations. The use of this mode for shear-wave velocity estimation has been investigated for many years (e.g., White, 1967; Chen, 1989; Ellefsen, 1990; Wang and Tang, 2003; Zhang et al., 2010).

Because borehole acoustic waves only penetrate a few wavelengths or borehole radii into the formation, it may be suspected that near-borehole effects could dominate sonic measurements significantly. In rocks where sonic speed is coupled to stress, this phenomenon is sometimes referred to as acoustoelasticity, variations in stress concentration around the borehole (i.e., hoop stresses) affect sonic velocities near the borehole. This effect has been clearly shown in cross-dipole logging (e.g., Sinha and Kostek, 1996; Winkler et al., 1998; Huang et al., 1999; Sinha et al., 2000). The potential impact of this stress effect on quadrupole sonic logging in a wellbore is the subject of this study. Stress concentration effects cannot be removed from the reservoir and measurement environment. From a modeling point of view, we may include them or leave them out of our analyses. Hence, modeling can quantify their effect, assess their importance, and revise log interpretations where needed.

Acoustoelastic coupling is considered to be due to stiffening of intergranular contacts of the hertzian type, or closure of microcracks or cracklike components of the pore space. For an acoustoelastic material, a nonhydrostatic state of stress causes orthorhombic symmetry; i.e., the elastic properties are those of an orthorhombic medium, with material axes aligned with the principal axes of stress. The constitutive model that we use to capture this stress-moduli interdependency is originally proposed by Mavko et al. (1995).

In our numerical study, we consider a formation with properties similar to Berea sandstone, for which experiments show that velocities increase with increasing confining stress (Sayers et al., 1990; Fang et al., 2013). We consider a circular vertical borehole, for which we characterize the perturbed stress and velocity field and simulate the quadrupole wave under heterogeneous stress and velocity conditions. We establish the dispersion characteristics and demonstrate how effects related to the near-borehole stress field can be accounted for and the true formation shear velocity can be estimated.

## FORMULATION

Propagation of acoustic waves in a borehole has been subject to extensive analytical works. A monograph around the topic has also been published (Tang and Cheng, 2004). Analytical methods work well in an unbounded linear elastic and isotropic homogeneous medium, but they have limitations when it comes to real formations. Real rocks are in many cases acoustoelastic, e.g., their velocities and stiffnesses increase with increasing confining stress (e.g., Nur and Simmons, 1969); hence, rocks are anisotropic and the formation is heterogeneous in the immediate proximity of the borehole. However, for reasons of modeling difficulty, the related impact of this effect on borehole acoustics is not fully implemented in most modeling methods.

We therefore derive a simulation model applicable for studying borehole acoustics in stress-sensitive formations that includes the elastic stress concentration around the borehole, and we assume orthorhombic material properties derived from a stress-dependent constitutive model (Mavko et al., 1995). We use a special purpose finite-element approach (Jørgensen, 1993), which has proven particularly useful in modeling stresses and strains in anisotropic elastic bodies of cylindrical geometry. The modeling approach was originally developed for materials science research, in which the study of composite materials has demanded numerical techniques for this particular class of problems (Jørgensen, 1994; Jørgensen et al., 1998).

### Constitutive model for stress-dependent elasticity

The constitutive model (Mavko et al., 1995) regards a class of rocks in which elastic closure of a random distribution of compliant cracks or cracklike voids is stiffening the material under increasing confining stress. The derived elasticity (compliance) tensor $Sijkl$ is isotropic (two independent parameters) under hydrostatic confining stress but orthorhombic (nine independent parameters), with material axes parallel to the directions of principal stresses, under a general triaxial loading.

Under the assumption that closure of microcracks, or cracklike components of the pore space, causes the elastic stiffening under increasing confining stress, the full stress-dependent compliance tensor $Sijkl(\sigma )$ is established. Here, $\sigma $ refers to the full stress tensor and the compliance tensor is derived for an arbitrary state of stress.

### 3D modeling of stresses and strains

Anisotropy and spatial variations in anisotropy make the elastostatic problem challenging. Some observations of symmetry are, however, very helpful. First, because we consider reservoirs under considerable burial, we assume that the vertical direction ($z$-direction) is one of the three orthogonal principal stress orientations. Second, we choose to align the global (*x*-, *y*-, *z*-) coordinates with the maximum and minimum horizontal stresses in the reservoir in all of our analyses, thus the global coordinate system ($x,y,z$) has axes parallel with the in situ stresses; hence, $\sigma xx,\sigma yy,\sigma zz$ are principal stresses.

We consider a vertical hole in the reservoir subject to in situ stresses ($\sigma xx,\sigma yy,\sigma zz$). We seek the stress distribution in the neighborhood of the hole, where a considerable stress concentration effect is present. The principal stresses and orthorhombic symmetry planes will have rotational symmetry about the vertical $z$-axis only as long as the two horizontal in situ principal stresses are equal. For all other cases, the elastic anisotropy contradicts the basis of axial symmetry. Still, however, the vertical planes of mirror symmetry, indicated in the horizontal ($r,\theta $-plane in Figure 1) are noted. The two planes of symmetry go through the borehole center and are parallel with each of the horizontal in situ principal stresses.

The benefit of the applied method is the reduction in computational effort compared to conventional 3D finite-element analyses, which is the result of a reduced number of degrees of freedom, and a reduced bandwidth of the system stiffness matrix relative to conventional 3D finite-element modeling. This has to do with the fast convergence of the series in equation 7, which in turn has to do with the relatively weak coupling of modes (equation 12).

### Sonic wave propagation along a borehole

All borehole-guided waves are eigenmodes and their shapes and associated frequencies can be derived as eigensolutions to a linear undamped eigenvalue problem. The eigenvalue problem is the wave equation formulated in the frequency domain excluding external forces. We now derive the eigenvalue problem and finally demonstrate a numerical formulation benchmark against known analytical solutions.

In equation 16, $Ks$ is the system solid stiffness matrix, $Kf$ is the system fluid stiffness matrix, $Ms$ is the system solid mass matrix, and $Mf$ is the system fluid mass matrix, $Q$ is the system coefficient of coupling, and $\omega $ is the angular frequency. The total number of degrees of freedom is $(1+N)Lf+(2+3N)Ls$, where $Ls$ and $Lf$ are the number of solid and fluid element nodes, respectively, and $N$ is the Fourier series truncation number.

### Test of numerical approach

The phase velocity is frequency dependent, and a plot of pairs ($f,Vphase$) reveals the dispersion characteristics (Figure 3). The simulated dispersion curve corresponds very well to the analytical counterpart; thus, the error related to discretization can be considered minimal.

In the present analysis, which is purely based on real number arithmetic, the cut-off frequency is identified as the frequency where transition to leaky modes occurs, i.e., where the pure shear deformation of horizontal polarization is excited, and the phase velocity approaches the true formation shear-wave velocity ($VS=G/\rho s$), which in this case is $2300\u2009\u2009m/s$. Mode shapes are shown for different frequencies, and the gradual change in mode shape toward pure shear near the cut-off frequency is displayed in Figure 4. Note that mode shapes are shown as dimensionless.

### Stress distribution around a vertical borehole

Simulation of the effect of the direction and magnitude of principal stresses on elastic moduli and dispersion curves provides the framework for understanding how velocities are perturbed near the borehole.

We use Berea sandstone for our study case. For this rock, a nonlinear stress-strain relationship of the type previously described in the section “Constitutive model for stress-dependent elasticity” applies, and we can use the $VP$ and $VS$ versus the $\sigma hyd$ relations published by Fang et al. (2013) (Figure 5). Because velocity and compliance are uniquely related, we can calculate the isotropic elastic compliances as a function of hydrostatic stress. Now, for any state of stress, which includes three principal stresses of different magnitude, we can derive the resulting full orthorhombic stiffness tensor as described in equations 1–6. The asymptotic compliance is approximated by the compliance at the highest hydrostatic stress point. This means the $VP$ and $VS$ velocities at 50 MPa (see Figure 5) are used to calculate the linear elastic compliances $Sijkl0,ISO$ and the trend curves (red and blue curves) are used to calculate compliances at lower confining stress $SijklISO$.

We consider a vertical hole subject to in situ stresses ($\sigma xx,\sigma yy,\sigma zz$). The stress distribution in the solid is established by imposing the in situ stresses on the outer circumference of the model $r=R$ and then solving the static problem. Stress-free conditions at the fluid-solid interface are applied, and the vertical displacements are fixed at $z=H$ and $z=0$. We truncate the Fourier series at $N=3$. The static boundary value problem is a nonlinear one due to the nonlinear stress-strain relationship, and the solution is derived using an incremental approach; i.e., we impose the in situ stresses in increments and solve the linearized problem at each loading step. At each loading step, the stiffness matrix is updated in accordance with the concurrent stress state. Importantly, the angular variation in stress, and therefore in the stiffness tensor, is taken into account (see Figure 1 and equation 13). At the last loading step, the distributed stiffness tensor in the model reflects the stress variation in the model as defined by the constitutive model for stress-dependent elasticity (equations 1–6).

### Velocity distribution and quadrupole wave dispersion

Because the orthorhombic elasticity tensor is uniquely related to the stress tensor and velocities are analytical functions of moduli, we may easily derive the perturbed velocity field around the borehole. In Figure 6, we compare the derived $VP1=C11/\rho s$ and $VP2=C22/\rho s$ to the principal stresses (Figure 6a and 6b).

Of particular interest is the shear velocity that comprises the upper limit of the quadrupole phase velocity. Unlike the homogeneous case, this limiting velocity is not known a priori when stresses and moduli vary in space. In the stress-sensitive medium, the slowest shear wave is one that is polarized in the direction of minimum stress and propagating along one of the two other orthogonal axes. The speed of the wave that is polarized in the direction of least horizontal compression ($i=2$) and propagating in the vertical direction ($i=3$) is $VS32=C44/\rho s$, whereas the speed of the wave propagating along the alternate horizontal material axis is $VS12=C66/\rho s$.

Far away from the wellbore, the stress conditions are homogeneous; e.g., $C44=C66<C55$ for the anisotropic case and $C44=C55=C66$ for the isotropic case. Near the wellbore, shear velocities are influenced by the stress concentration (see Figure 8).

The spatial distributions of $VS12$ and $VS32$ for the anisotropic (a and b) and isotropic (c and d) in situ stress conditions are shown in Figure 8. For the isotropically stressed formation, the S-wave velocity attains its minimum at the borehole wall at $2300\u2009\u2009m/s$, which is to be compared to $2355\u2009\u2009m/s$ farther away from the borehole.

As expected, the variation is much greater for the anisotropic case, where the range is as wide as $2100\u20132400\u2009\u2009m/s$ for either of $VS12$ and $VS32$. For comparison, the horizontal P-wave velocities computed for isotropic in situ stress conditions vary between 3540 and 3750 m/s (radial and circumferential directions, respectively, see Figure 7c and 7d) and the horizontal P-wave velocities computed for anisotropic in situ stress conditions vary between 3500 and 3900 m/s. We now direct our attention to how the quadrupole phase velocity and dispersion characteristics are affected by the perturbed velocity field. To this end, we again focus on the eigenvalue problem (equation 16), and again we solve it repeatedly for different wavenumbers $kz$. That is, we alter the vertical dimension of the model step by step and solve the boundary value problem (equation 16) using the boundary conditions (equation 37) and the symmetry condition (equation 38). Clearly, the solid stiffness matrix is evaluated by solving exactly the static problem described above, and the remaining matrices in the formulation of equation 16 are evaluated as described in the section “Sonic wave propagation along a borehole.” Each eigenpair comprises a mode and a frequency, and we focus consistently on the quadrupole eigenmode (see Figure 2a). In a strict sense, the quadrupole eigenmode is not characterized by the acoustic profile $cos(2\theta )$ due to coupling related to anisotropy. The term “quadrupole eigenmode” refers, in the anisotropic case, to the eigenmode dominated by the azimuthal acoustic profile $cos(2\theta )$. Numerically, the quadrupole eigenpair (eigenmode and frequency) is again found by the iterative algorithm known as inverse iteration, this time with a shift. One advantage provided with this method is the convergence to the eigenvalue nearest to a predefined estimate (the shift value). An approximation to the quadrupole eigenvalue is easily obtained by truncating the series, given in equation 7, at $N=1$, and eliminating the degrees of freedom associated to the axial symmetric profile $n=0$ from the eigenvalue problem (equation 16). In this way, a reduced eigenvalue problem is obtained, i.e., one where coupling is eliminated; however, the lowest eigenvalue for each wavenumber $kz$ is an approximation to the quadrupole squared eigenfrequency. To each quadrupole eigenpair, we associate the mode’s phase velocity $Vphase=\omega /kz$ and frequency $f=\omega /2\pi $ and plot the dispersion curve ($f,Vphase$). These curves are plotted in Figure 9 for the isotropic as well as for the anisotropic cases. Also in Figure 9, the results of simulations are shown in which the stress concentration effect on the borehole hoop stress has been removed from the analysis. For the latter cases (dotted lines), the velocity fields are homogeneous and reflect conditions far away from the borehole. The curves converge toward the true formation shear-wave speed $VS32=VS12$ as expected. The curves influenced by near-borehole stresses are consistently below their undisturbed counterparts, except at the cut-off frequency where they coincide. Intuitively, this is to be expected because in either cases, the shear mode is excited in the entire solid phase and the speed of propagation is $VS32=VS12$. However, it should be noted that shear velocity estimates in field measurements are based on extrapolations of quadrupole dispersion curves to the cut-off frequency because the mode is generally excited above the cut-off frequency only. To this end, we note that for the anisotropic case, the influence of the borehole hoop stress concentration causes the phase velocity trend to be around 1% less than without the stress influence. This would imply that the quadrupole dispersion is essentially insensitive to near-wellbore effects for the cases considered, even when the spatial variation of shear velocities is quite significant near the borehole. This is consistent with approximations derived from perturbation analyses (Sinha et al., 2003) and recent experimental results (Hsu et al., 2011).

Information about anisotropy can be extracted from analysis of quadrupole modes. The quadrupole mode, which is associated with azimuthal variation $cos(2n\theta )$, couples to modes with azimuthal variation $cos(4n\theta )$, $cos(6n\theta )$…, as well as to axisymmetrical modes. The strength of coupling reflects the degree of anisotropy. Also, the azimuthal directions ($\theta ,\theta +\pi /2$) of maximum amplitude coincide with the directions of maximum and minimum in situ horizontal stress.

## CONCLUSION

We have derived an integrated numerical finite-element approach, which includes couplings between nonlinear elastic deformation, near-borehole stress redistribution, and borehole acoustics. The coupling between stress and elasticity means that spatial variations in elastic anisotropy apply near a borehole, where stress concentration effects are present. We demonstrate how the method is constructed to include heterogeneous anisotropy effects for the case of a vertical wellbore in a stress regime where two in situ principal stress directions are horizontal and one is vertical. We have used the method to study the dispersion characteristics of quadrupole waves traveling along a borehole. Our analysis suggests that the dispersion characteristics of quadrupole waves are essentially insensitive to the spatial variation of elastic anisotropy immediately around the borehole. Hence, extrapolation of a quadrupole dispersion curve may provide estimates of the true formation shear velocity, independently of the spatial variation of this parameter in the close proximity of the borehole. This finding, which is in line with experimental results and approximations derived from perturbation analysis, is here confirmed in a full 3D analysis.

The presented method is based on a ring-type finite element, thus a radial/axial plane is taken as a plane of reference. The field variation out of this plane, i.e., the displacement field in the rock and the pressure field in the fluid, is restricted to a subspace formed by a truncated Fourier series. Full 3D analysis, at a computational effort that is remarkably reduced compared to conventional 3D finite-element modeling, can thus be accomplished for a small but very important class of problems by the proposed method. It is noteworthy that this class of problems includes the possibility of taking the effect of the bottom of the hole into account in the modeling. Hence, the method has potential for investigating applications for imaging ahead of the bit during drilling.

## ACKNOWLEDGMENTS

The authors wish to thank Maersk Oil and Gas AS for sponsoring this work and the Earth Resources Laboratory (ERL) at MIT for hosting it. The authors also wish to thank the students, staff, and faculty at MIT-ERL for many fruitful discussions.