The characterization of fractures is essential to increase the production of hydrocarbon and geothermal resources. In this study, we investigate the effect of the fluid flow in the longitudinal fracture along a borehole on the dispersion and attenuation behavior of Stoneley waves using numerical experiments. In general, incorporating a fracture with the aperture of several tens to hundreds of micrometers into 3D seismic modeling is challenging with high calculational costs. We develop a novel numerical scheme that includes a 2D fluid flow simulation embedded into a 3D wave propagation modeling to address this problem. We devise an approach for capturing the effects of the fluid flow in the fractures of arbitrary aperture widths, which could be much thinner than the grid spacing of a 3D wave propagation simulation. A comparison of the numerical results from our scheme with the analytical solution indicates good agreement, supporting the method’s validity. This developed scheme is applied to the coupled simulation between the Stoneley wave propagation along the borehole axis and induced fluid flow inside a longitudinal fracture with different fracture apertures, fluid viscosity, and dynamic hydraulic conductivity. The modified matrix pencil algorithm applied to the recorded waveforms estimates the dispersion and attenuation of the Stoneley mode. The numerical results reveal the effect of the fracture aperture and fluid viscosity on the dispersion and attenuation behavior of the Stoneley waves. Based on the results, we demonstrate our scheme is an innovative method for estimating the aperture of the fracture and viscosity of the fluid by analyzing the dispersion and attenuation properties of the Stoneley waves.

Fractures are subsurface conduits that allow fluid to pass through them more efficiently than through porous formation matrices (Carey et al., 2020). In the production of hydrocarbon and geothermal resources, hydraulic fracturing techniques are applied to create fractures to improve productivity (Bakku et al., 2013). The fractures induced by hydraulic fracturing or drilling propagate radially from the longitudinal openings formed on the opposite sides of the borehole wall in the direction of the maximum principal stress. Fracture geometry and permeability affect the propagation of borehole mode waves along the borehole walls, such as the flexural and Stoneley waves. Zheng et al. (2009) study the influence of a drilling-induced tensile fracture on cross-dipole measurements. Kayama et al. (2021a, 2021b, 2021c) investigate the effect of the depth of the longitudinal fracture on the flexural dispersion curve. Lei and Sinha (2013) study fracture-induced anisotropy and the reflectivity of the Stoneley and flexural waves at a fluid-filled fracture in the borehole. Tang and Cheng (1989) investigate dynamic fluid flow in the planar parallel fracture. Tang et al. (1991) investigate the Stoneley dispersion and attenuation in the fluid-filled fractures in a borehole exhibiting multiple fracture openings using an analytical approach. In the approach that Tang et al. (1991) adopt, the viscous skin depth is small compared with the fracture opening. The interaction between the Stoneley wave propagating along a borehole axis and fluid flow inside the longitudinal fracture (we call this interaction a “borehole-fracture interaction” in the following) is a direct consequence of the viscous resistance of the fluid at the two fracture faces (Gehne and Benson, 2019). Tang et al. (1991) point out the effects of a borehole-fracture interaction increasing at low frequencies (less than approximately 0.5 kHz) in the propagation of the Stoneley waves. Considering the possible variation of a fracture aperture (micrometer order) or signal frequencies of the Stoneley waves (kilohertz order), we think a quantitative dispersion analysis of the Stoneley waves needs to be conducted in a nonnegligible viscous skin depth case.

The influence of fluid viscosity in a fracture radially intersecting a borehole has been studied by Bakku et al. (2013). They use numerical analysis of dispersion equations developed by Tang (1990) for the viscous skin depth comparable to or larger than the fracture opening. Because the target aperture has been in the order of 1 mm or broader (Tang, 1990), the viscous effects are not included in the case of longitudinal tensile fracture because the dynamic conductivity is strongly dominated by the propagation of the acoustic wave in the fracture (Tang and Cheng, 1989). One of the schemes to estimate the influence of a thin fracture filled with a viscous fluid is to accommodate a possible interaction with a numerical simulation. The boundary conditions that define the interaction between the formation and viscous fluid need accordingly to include the viscous resistance of the fluid because the viscous shear effects attenuate the energy of the propagating waves (Schoenberg et al., 1987). If we discretize the whole computational domain (borehole, fluid-filled fracture, and surrounding formation) using a regular lattice, the degrees of freedom reach 1012–1013 order or more when the fracture aperture is less than 1 mm. In this case, the required computational memory is more than several tens or hundreds of terabytes at a rough estimate. Therefore, for numerical simulation, mesh or particle generation to define the solid-fluid boundary of a thin fracture requires substantial computational resources. Because the computational model for 3D seismic wave propagation accepts coarser nodal spacing, separation of 3D seismic wave propagation and 2D fluid flow in the longitudinal fracture would be a helpful strategy in some way. Although the effective medium theory can incorporate the effects of small fractures, which are negligible compared with the wavelength of the propagating seismic motion, into macroscopic physical properties, it is difficult to consider the effects of fluid flow in a fracture on propagating seismic waves. For estimating the effects of fluid flow on the Stoneley wave propagation, it is necessary to introduce a new method to capture the effect of fluid flow without sacrificing computational efficiency.

This paper, referring to the dynamic fracture flow theory (Tang and Cheng, 1989), performs the coupling calculation of the 2D fracture flow and 3D sonic logging. The effect of the aperture of the longitudinal fracture and viscosity of fluid on Stoneley wave dispersion and attenuation is investigated. Because there is a large gap of several orders of magnitude between the diameter of the borehole and the opening width of the small-aperture fracture, sonic logging simulation in this environment becomes a multiscale problem. We develop a novel numerical scheme for simulating the solid-fluid coupling effect. The developed scheme is verified by comparing it with the analytical dispersion curve of the fracture mode by Ferrazzini and Aki (1987). Then, a 3D borehole model including longitudinal fractures is created, and a monopole sonic logging simulation is performed with different fracture apertures, fluid viscosity, and dynamic hydraulic conductivity of various representative frequencies (we will define “representative frequency” in the following section). In addition, we show a possibility for estimating the fracture aperture and fluid viscosity by analyzing the dispersion and attenuation properties of the Stoneley waves from the numerical results.

This section explains the proposed numerical scheme for simulating seismic wave propagation, fluid flow, and their coupling effect.

3D seismic wave propagation

The Hamiltonian particle method (HPM) (Takekawa et al., 2014a) simulates 3D seismic wave propagation in the present study. The HPM can estimate seismic wavefields accurately in numerical models with complex boundaries (Takekawa et al., 2014b). In the present study, we need to discretize the model around the borehole wall with fine nodes because higher frequencies are affected by small steps at the solid-fluid boundary (Takekawa et al., 2014b). Furthermore, low-frequency components induce the displacement field in the broader area (i.e., deeper from the borehole wall) more than high-frequency components. Therefore, the absorbing boundary should become thick enough away from the borehole to avoid any incorrect attenuation of low frequencies. Substantial computational costs are required if we use a regular latticed node for discretization. Therefore, we introduce a nonuniform distribution of particles into the HPM to assure accuracy in low- and high-frequency components without sacrificing computational efficiency. Details of the method are shown in Appendix  A. Although the flexural wave propagation by the HPM was validated by comparing the analytical dispersion curve and the waveforms by the discrete wavenumber integration method in Kayama et al. (2021a), additional validation for the Stoneley wave propagation by the HPM with nonuniform particle distribution is required. The validation result also is shown in Appendix  A. The modified matrix pencil algorithm (Ekstrom, 1995) then estimates the dispersion and attenuation from the obtained shot gathers.

2D fluid flow

Fluid flow in the longitudinal fracture is modeled by 2D fluid flow in a parallel plate. The dynamic hydraulic conductivity in the fracture is governed by two effects, i.e., diffusion and propagation, which are determined by the signal frequency, fracture aperture, and fluid viscosity. The basic theory is based on Tang and Cheng (1989). In the following, we explain the theory of 2D fluid flow.

The dynamic hydraulic conductivity C (Tang and Cheng, 1989) following Darcy’s law is
q=Cp,
(1)
C=iωLk2αf2ρf,
(2)
where q is the unit length flow rate, p is the pressure gradient (averaged within fracture aperture), ω is the angular frequency, L is the fracture aperture, k is the wavenumber, αf is the fluid acoustic velocity, and ρf is the fluid density. The wavenumber k is the solution of the following characteristic equation representing a fracture mode that propagates along the fracture between two parallel plates:
k2tan(f¯L2)+ff¯tan(fL2)=0,
(3)
f2=ω2αf243iωνk2,f¯2=iωνk2,
(4)
where ν is the fluid kinematic viscosity. As indicated by Tang and Cheng (1989), the dynamic hydraulic conductivity becomes diffusion dominant (C=L3/12μ: cubic law) in the low-frequency/low-aperture/high-viscosity limit and propagation dominant (C=iL/ωρf) in the high-frequency/large aperture/low-viscosity limit. From equations 2 and 3, the dynamic hydraulic conductivity depends on the frequency of propagating seismic waves. Therefore, the dynamic hydraulic conductivity in the numerical simulation should be changed by the signal frequency. However, in the time-domain method, the dynamic hydraulic conductivity should be set at a certain frequency, whereas the source wavelet generally includes wide frequencies. Here, we define representative frequency as a frequency f0(=ω0/2π) which determines the dynamic hydraulic conductivity in equation 2. The numerical simulation is conducted with the designated representative frequency, and then the dispersion and attenuation properties are calculated using recorded waveforms. Subsequently, we extract the slowness and attenuation coefficient from the dispersion and attenuation curves at the representative frequency. The dynamic hydraulic conductivity is decomposed into real and imaginary parts:
C=iω0Lk2αf2ρf=kcrL312μ+ikciLω0ρf,
(5)
where μ(=ρfν) is the viscosity of the fluid. Here, the coefficients kcr and kci are constants that describe the effects of diffusion and propagation, respectively. The coefficients are determined for any given aperture and fluid viscosity by fixing the representative frequency f0. The following section will show the relationship between the coefficients and parameters (aperture, fluid viscosity, and representative frequency). The average flow velocity v in the fracture is
v=qL=CLp=kcrL212μpikci1ωρfp.
(6)
Transforming this equation into the time domain leads to the following equations:
dζdt=kci1ρfp,
(7)
v=kcrL212μp+ζ,
(8)
where t is the time and ζ is a memory variable that describes the flow velocity due to propagation. We simulate fluid flow in the fracture by discretizing these equations. The derivation of equations 7 and 8, discretization, and their validation are shown in Appendix  B.

Coupling of 3D seismic wave propagation and 2D fluid flow

This section explains the coupling strategy between 3D seismic wave propagation and 2D fluid flow simulations. Figure 1 shows the spatial distribution of variables in the vicinity of the longitudinal fracture. This calculation uses two types of particles, i.e., stress and velocity particles. Stress and pressure are defined at the stress particles, whereas velocity is defined at the velocity particles. Physical properties are defined at the stress particles. Density must be defined at the location of velocity particles and should be interpolated by the surrounding stress particles (Saenger et al., 2000). The Lamé constants are used only at the stress particles; therefore, such interpolation is unnecessary. Numerical calculations for fluid flow in the fracture are called a “2D system,” whereas those for seismic wave propagation are called a “3D system.” Additional velocity particles for the 2D system are added between the stress particles located on the fracture plane. The stress particles are shared by 2D and 3D systems. The results of the two systems are coupled at every discrete time step in the simulation. Because 2D and 3D systems are coupled after simulating them individually, arbitrary apertures of the fracture can be implemented independently with respect to the particle spacing in the 3D system.

Figure 2 shows a time flow chart of the calculation process. The pressure in the fracture is updated in the 3D and 2D systems. However, velocities in the 2D and 3D systems are updated individually. In terms of the stability condition, the diffusion calculation in the 2D system has the strictest stability condition. If we set the time step to coincide with the stability condition of the diffusion calculation, the total calculation cost increases drastically due to the seismic wave propagation in a 3D system. Therefore, two types of time steps are used in this study, i.e., one for the diffusion calculation in a 2D system and the other for the propagation calculation in a 2D system and seismic wave propagation in a 3D system. Because the time step for the diffusion calculation is finer than the other, several iterations are required during one time step in the propagation calculation in the 2D system and the seismic wave propagation in the 3D system. This time updating strategy will be validated in the next section.

For the seismic wave propagation in a 3D system, the physical properties (density and bulk modulus) of particles at the fracture location require special treatment. We set their physical properties in the following manner. Figure 3 shows a schematic figure of a 3D elementary volume, including the longitudinal fracture. The elementary volume consists of the fluid and solid parts, and the fracture passes through the stress particle. The volume fraction of the fluid part is defined as f. The volume of the fluid part Vf and solid part Vs can be represented by f and the volume of composite volume Vc as follows:
Vs=(1f)Vc,
(9)
Vf=fVc,
(10)
where subscripts “f,” “s,” and “c” represent the fluid, solid, and composite parts, respectively. The volume change rate dΔVi can be written as follows:
dΔVi=deiVi,(i=s,f,c),
(11)
where dei is the volumetric strain rate. The density of each part is calculated based on the mass conservation law:
mc=ms+mf,
(12)
ρcVc=ρsVs+ρfVf,
(13)
ρc=(1f)ρs+fρf,
(14)
where mi and ρi are the mass and the density, respectively. Next, we calculate the bulk modulus Kc that is averaged over the elementary volume as in effective-medium theories (e.g., Yue and Yue, 2017). The pressure change dpi is equal in the fracture and formation parts:
dpc=dpf=dps.
(15)
At this time, the sum of the volume change rates dΔVi of the solid and the fluid parts are the same as the volume change rate of the composite part:
dΔVc=dΔVs+dΔVf,
(16)
decVc=desVs+defVf,
(17)
dpcKcVc=dpcKs(1f)Vc+dpcKffVc.
(18)
That is,
1Kc=1fKs+fKf.
(19)
This equation agrees with the Reuss average for the effective fluid bulk modulus of a partially saturated porous medium (Krief et al., 1990) when the dry bulk modulus is negligible compared with the matrix bulk modulus using the Biot-Gassmann theory (Mavko and Mukerji, 1998).
We assume that the volumetric strain def0 (or dpf0) generated by the fluid flow calculation in the 2D system is instantly mixed over the composite volume as shown in Figure 4. This assumption is valid as long as the propagating wavelength is sufficiently longer than the nodal spacing in the opening direction. This assumption always holds in the following numerical experiments. Then, the pressure change in the elementary volume after mixing dpc can be calculated. The volume change rate dΔVi before and after the mixing is the same as each other:
dΔVf0=dΔVc,
(20)
def0Vf=decVc,
(21)
def0fVc=decVc,
(22)
def0f=dec.
(23)
At this time, the pressure change dpc after mixing is
dpc=Kcdec=fKcdef0.
(24)
Therefore, the coefficient Kfc for calculating the pressure change dpc is as follows:
Kfc=dpcdef0=fKc.
(25)
Pressure change before the instant mixing dpf0 is
dpf0=Kfdef0.
(26)
Therefore, the ratio fc of the pressure change of the fluid part before and after the instant mixing is
fc=dpcdpf0=Kfcdef0Kfdef0=KfcKf.
(27)
However, in the whole fluid grid,
Kfc=dpfdef0=fKf,fc=f.
(28)

Validation of the proposed scheme

For the validation of the coupling theory, we use a simple 3D numerical model with a fracture, as shown in Figure 5. Particle spacing is Δx=Δz=6.25  mm,Δy=1  mm around the fracture. The spacing extends with increasing the distance from the fracture. The z-direction is extended infinitely by a periodic boundary condition. A 3 kHz dipole source is located on the fracture (3D velocity particle of both sides) and pressure in the fracture is recorded by a receiver array. In this case, the symmetric mode (sausage mode) is induced inside the fracture. We chose this source type for validation because the symmetric mode would be induced in the longitudinal fracture by the Stoneley wave propagation in sonic logging. The physical properties of formation and fracture are shown in Table 1. Conductivity in the fracture is set to a propagation-dominant condition (kcr=0,kci=1). The fracture aperture L=50  μm is much smaller than the grid spacing Δy (1/20). In this case, we can use the analytical dispersion equation of the fracture mode (Ferrazzini and Aki, 1987) as follows:

cot(Lω2υυ2αf21)=ρVS4υ2αf21ρfυ41υ2VP2[(2υ2VS2)241υ2VS21υ2VP2],
(29)
where υ is the phase velocity of the fracture mode, ρ is the density of the formation, VP is the P-wave velocity of the formation, and VS is the S-wave velocity of the formation. Tang and Cheng (1988) show that the experimental fracture mode is consistent with this equation.

Figure 6 shows the snapshot of the y-directional pressure, i.e., a shot gather recorded by the receiver array and the dispersion curves for the gather. The dispersion curve (the red line) is estimated using the modified matrix pencil algorithm, and the spectral semblance of the numerical waveform also is plotted as the background counter. Numerical dispersion is in perfect agreement with the analytical dispersion curve. This result shows that the proposed coupling method can accurately reproduce the velocity dispersion of the fracture trap mode in that the aperture is much smaller than the grid size.

As described in the previous section, we test the effect of the difference between the two time steps, one for the diffusion calculation in the 2D system and the other for the propagation calculation in the 2D system and seismic wave propagation in the 3D system. Three patterns of time steps shown in the following are examined:

  1. ζ is updated in dt3D (not dependent on mixing timing),

  2. ζ is updated in dt2D, instant mixing in dt2D,

  3. ζ is updated in dt2D, instant mixing in dt3D,

where dt3D=2.5×107  s and dt2D=2.0×109  s=dt3D/125.

Here, dt2D is a time step for the diffusion calculation in the 2D system and is set to 125 times finer than dt3D. The test model is the simple fracture model shown in Figure 5. Figure 7 shows the attenuation and dispersion curves with each pattern. Although the dispersion curves agree with each other, the attenuation property shows significant differences. Because there is no attenuation in this fracture mode, the time updating strategy of pattern 1 is appropriate. In other words, patterns 2 and 3 generate false attenuation. Although the reason why patterns 2 and 3 generate artificial attenuation has not been clarified in this study, Figure 7 clearly shows that pattern 1 can eliminate artificial attenuation. Therefore, in the following section, we choose the time steps of pattern 1, i.e., 2.0×109  s for the diffusion calculation in the 2D system and 2.5×107  s for the propagation calculation in the 2D system and seismic wave propagation in the 3D system.

The preceding results are obtained under the propagation-dominant condition. Next, simulations of the fracture mode propagation are performed at various representative frequencies with a viscosity of μ=4  mPa·s. Numerical experiments are carried out with the following three time-setting patterns:

  1. dt2D=t3D (mixing timing also is the same),

  2. dt2D=2.0×109  s=dt3D/125, dp is mixed immediately in dt2D,

  3. dt2D=2.0×109  s=dt3D/125, dp is mixed in dt3D,

where dt3D=2.5×107  s, and ζ is updated in dt3D (pattern 1).

The results from the preceding conditions are compared with the analytical dispersion and attenuation curves obtained from the theoretical dispersion equation for the fracture mode involving viscous fluid. The derivation is shown in Appendix  C. Figure 8 shows the comparison results between the numerical and analytical solutions. For the numerical solution, the results for the time-setting patterns A–C are shown. The results obtained from all calculation conditions show good agreement for the dispersion curves. As for the attenuation curve, the time-setting pattern A shows almost perfect agreement, whereas patterns B and C show some excess attenuation. The results of pattern A indicate that the proposed coupling method can accurately reproduce the dispersion characteristics of the fundamental fracture mode with viscous fluid when the same time steps are adopted for the 2D and 3D systems. However, it is difficult to use the same time steps for 2D and 3D systems in terms of the computational costs due to the strict stability condition of the diffusion equation for high permeability (L=100  μm, μ=4  mPa·s) with fine particle spacing (1 mm) near the borehole. Therefore, pattern B is used in the following numerical experiments because the error in the attenuation is small enough (less than 3%), i.e., we choose 2.0×109  s for the diffusion calculation and mixing timing of the 2D system and 2.5×107  s for the propagation calculation of the 2D system and seismic wave propagation of the 3D system.

We conduct the sonic logging simulation using the proposed coupling method. Figure 9 shows the numerical model geometry for the fluid-filled borehole with a longitudinal fracture. The borehole radius is set to 10 cm. Particle spacing is Δx=Δy=1  mm, Δz=6.25  mm around the borehole. The spacing in the x- and y-directions is expanded gradually with increasing the distance from the borehole. The surrounding formation is fast formation, whose physical properties are shown in Table 1. The 3 kHz monopole source whose waveform is the Ricker wavelet is located in the center of the borehole. Forty receivers are placed at 5 cm intervals from 1 m offset above the source. The waveform recording time is 15 ms.

We conduct numerical experiments with a different aperture of fracture and viscosity of the fluid. Table 2 shows a combination of them. Both parameters have three patterns, and the dynamic conductivity of the fracture, which is determined by the aperture of fracture and viscosity of the fluid, also has three patterns by characteristic equation 3. Figure 10 shows the variations in the contribution of diffusion and propagation with representative frequency for each numerical experimental condition.

Figure 11 shows an example of the simulation result in the propagation-dominant condition (L=50  μm, kcr=0, kci=1). We can see that the fracture mode propagates toward the outside of the model. Because the fracture mode does not return to the borehole, it causes an energy loss of the borehole Stoneley mode. The flux of dissipated energy flows out not only from the fracture fluid but also from the fracture wall to the formation part.

At first, we show the effect of the fracture and fluid flow on the propagation behavior of the Stoneley wave. Figure 12 shows snapshots of the pressure field in different models, i.e., the borehole model without a fracture, the borehole model with a fracture but the fluid flow is not considered, and the borehole model with a fracture and the fluid flow in it. We can observe the leaking wave along the fracture in the model with the fracture and fluid flow, whereas the wave energy is trapped in the other models. Figure 13 shows the dispersion and attenuation curves for each model. The dispersion property is almost the same for each model, although a slight difference can be observed in the low-frequency region. However, the attenuation property shows a clear difference with or without the fluid flow. As shown in Figure 13, the fluid flow in the fracture causes the attenuation of the Stoneley wave because the lower frequency induces much more fluid flow inside the longitudinal fracture. This trend is similar to the previous study (Tang et al., 1991).

Next, we simulate the wave propagation with a different aperture of the fracture. The snapshots of the pressure field are shown in Figure 14. The fracture mode is induced by the wave propagation along the borehole in each model. It can be observed that the phase velocity of the fracture mode increases with increasing the fracture aperture. Because this result agrees well with the previous research (Tang et al., 1991), it can be confirmed that the proposed method can reproduce a reasonable result with different apertures. Although the investigation of Tang et al. (1991) is limited to millimeter- and centimeter-order apertures, the result extends its limitation to smaller apertures (i.e., micrometer order). Figure 15 shows the dispersion curves with different apertures. At the low-frequency region, the slowness increases (phase velocity decreases) with increasing the aperture. However, the dispersion curves with small apertures have almost no difference from those without fractures. The dispersion curve with the aperture of 1000  μm begins to move away from that without fracture at approximately 1 kHz, whereas the slowness with the aperture of 100  μm increases rapidly at approximately 0.2 kHz. This result is consistent with the study of Tang et al. (1991), whose theory agrees well with a laboratory experiment, i.e., showing the effectiveness of the proposed method. Because the change in the aperture only affects the dispersion properties at low frequencies, the slowness change at the low-frequency region would help us detect fractures with a large aperture. However, it is challenging to estimate small apertures using the dispersion property of the Stoneley wave.

Next, we investigate the attenuation property of the Stoneley wave. Figure 16a shows the attenuation curves with different apertures. All curves show apparent differences from each other in a wide frequency range due to the fracture fundamental mode at low frequencies (Ferrazzini and Aki, 1987; Tang and Cheng, 1988) discussed in Tang et al. (1991). Figure 16b shows the attenuation as a function of the aperture with different frequencies. The lower frequency components have higher sensitivity than the higher frequency components. Each attenuation curve with respect to the aperture can fit the power function. This result indicates that the aperture can be estimated by analyzing the attenuation property depending on the frequency. It should be noted that the kink of the attenuation curves at approximately 7 kHz caused by the nonuniform nodal distribution around the borehole is observed because the change of the nodal distribution pattern shifted the kink. The change of the nodal spacing may generate tiny reflection and refraction (they are invisible in snapshots, e.g., Figure 17). They could interact with the Stoneley wave propagation, and then the kink of the attenuation curve would be observed. In our model, this interaction has a significant effect, especially in a lower frequency, as shown in Figure 18b. Because the attenuation coefficient is picked at the representative frequency, we think that the effect of this false attenuation on the result can be negligible.

Next, we simulate the Stoneley wave propagation with different representative frequencies, i.e., the effects of propagation and diffusion are not zero. Parameters of the fracture and fluid are fixed to L=50  μm, μ=4  mPa·s. Figure 17 shows snapshots of the pressure field with the representative frequency ranging from 0.5 to 8 kHz in addition to propagation- and diffusion-dominant conditions. In low representative frequencies, the fracture mode shows strong attenuation, whereas its behavior approaches the propagation-dominant condition (Figure 11) with increasing representative frequency. Therefore, the diffusion effect in the longitudinal fracture has consequences for attenuating the Stoneley wave in the lower frequency region because the diffusion effect increases with decreasing the representative frequency. As mentioned previously, the dynamic hydraulic conductivity is fixed after determining the representative frequency. At first, we need to calculate dispersion and attenuation curves with various representative frequencies, as shown in Figure 18. Each curve is drawn by the simulation results with different representative frequencies. Then, we pick the slowness and attenuation coefficient at the representative frequency from the corresponding curves (the filled triangles in Figure 18). Each triangle is spliced, and a black broken line shows the dispersion and attenuation curves. Figure 18 is obtained with parameters of L=50  μm, μ=4  mPa·s. The peak of the attenuation curve tends to shift to the low-frequency side with increasing the representative frequency. Therefore, the magnitude of the attenuation at a higher frequency (in this case, higher than 3–4 kHz) is reversed for different representative frequencies. Next, we change these parameters to produce the dispersion and attenuation curves in the same manner.

Figure 19 shows the dispersion curves with different apertures and viscosities. Although the slowness at the low-frequency region slightly changes, little difference between each curve can be seen, which indicates that the dispersion property of the Stoneley wave is not sensitive to the fracture aperture. Although we consider the effect of the viscosity of fluid on the propagation behavior of the Stoneley wave, the trend of the dispersion curve is almost the same as in the previous studies. It is difficult to estimate the fracture aperture and fluid viscosity by analyzing the dispersion property of the Stoneley wave. Figure 20 shows the attenuation curves with different apertures and viscosities. In contrast to the dispersion property shown in Figure 19, apparent differences, especially in the low-frequency region, can be distinguished. At the high-frequency region, the fluid’s viscosity has a smaller effect on attenuation. The tendency in the attenuation can be confirmed by comparing the attenuation curves with the same colors in Figure 20. Therefore, the attenuation at the high-frequency region is mainly governed by the fracture aperture. The high frequency’s attenuation behavior could help us to estimate the aperture of the fracture free of influence from the fluid viscosity.

We have shown the possibility of estimating the aperture of the fracture and the fluid’s viscosity. However, all field data include noise to some extent. Next, we investigate the effect of noise on the estimation of the fracture apertures and fluid viscosities. Random noise is added to the synthetic data with different signal-to-noise ratios (S/Ns) (2, 5, and 10). We estimate the dispersion and attenuation curves using noisy data, as shown in Figure 21. This example is calculated with L=50  μm, μ=4  mPa·s, and an S/N of five. Although the curves are disturbed by the noise, their trend can be recognized at lower frequency ranges. However, it is difficult for higher frequencies to recognize the slowness and the attenuation coefficient. Similar results can be obtained for other calculation conditions. Therefore, we extract the slowness and the attenuation coefficients at the representative frequencies up to 6 kHz. Figure 22 shows the dispersion and attenuation curves with various S/N. Compared with Figures 19 and 20, the estimation of the dispersion curves is relatively robust. However, the estimation results of the attenuation are largely affected by the noise. This result indicates that the noisy environment may cause a misinterpretation of the aperture and fluid viscosity.

In the preceding numerical experiments, we find the following remarks.

  • Low-frequency slowness rises in large aperture cases.

  • The attenuation property is affected only by the fracture aperture in propagation-dominant conditions.

  • High-frequency attenuation is almost unaffected by the viscosity of the fluid.

  • The viscosity of the fluid suppresses low-frequency attenuation (following the characteristic equation).

  • Fractures with a small aperture and diffusion-dominant conditions induce less attenuation.

Based on the preceding findings, we propose a new scheme for estimating the aperture and viscosity using the dispersion and attenuation properties of the Stoneley wave. We assume that the attenuation curve is a function of three parameters (frequency, aperture, and fluid viscosity), i.e., αt(ω,L0,μ). If attenuation curves of the Stoneley wave can be calculated analytically for the given parameters, the aperture and viscosity are estimated by matching the observed attenuation curves with the analytical ones. However, the lack of symmetry makes it difficult to use an analytical approach for the longitudinal fracture case. We consider a linear independent model of the propagation- and diffusion-dominant attenuation curve based on the characteristic equation’s dynamic permeability coefficient:
αt(ω,L,μ)kci(ω,L,μ)×αp(ω,L)+kcr(ω,L,μ)×αd(ω,L,μ),
(30)
where αt(ω,L,μ), αp(ω,L), and αd(ω,L,μ) are the attenuation curves in transition frequency, the propagation-dominant conditions, and the diffusion-dominant conditions, respectively. Here, kcr(ω,L,μ) and kci(ω,L,μ) are the coefficients of the real part (diffusion) and the imaginary part (propagation) of conductivity, respectively.

To test the validity of the preceding equation, we calculate the attenuation curves with three different dynamic conductivities. We change the fluid viscosity (1, 4, and 16 mPa ⋅ s) with a constant aperture of 50  μm. We calculate the attenuation curves under propagation- and diffusion-dominant conditions with these apertures and viscosities. They are shown in Figure 23 as broken and dotted lines. For each frequency, the estimated attenuation curves are calculated in equation 30 and shown in Figure 23 as blue lines. The red lines are calculated attenuation curves that come from Figure 20. Although some errors can be seen around the transition zone from a diffusion-dominant condition to a propagation-dominant condition, the calculated (the red) and estimated (the blue) attenuation curves have an excellent agreement with each other. This result indicates that equation 30 is valid for estimating the attenuation behavior of the Stoneley wave.

Because the attenuation curve is a function of the aperture and fluid viscosity, it is difficult to estimate both individually. Therefore, we first estimate the aperture of the fracture at the high-frequency region because the attenuation behavior is almost unaffected by fluid viscosity at the high-frequency region. After determining the aperture of the fracture, the remaining parameter (fluid viscosity) can be estimated by fitting the observed attenuation curve to the estimated one as follows:
αt(ω,μ)kci(ω,μ)×αp(ω)+kcr(ω,μ)×αd(ω,μ).
(31)

The solution of the characteristic equation 3 and the attenuation curve in the diffusion-dominant conditions are monotonous functions with respect to viscosity, as shown in Figure 24.

This study proposed a coupling method between 3D elastic wave propagation and the associated 2D fluid flow in the monopole sonic logging situation. In the new coupling strategy, the interaction between the Stoneley wave propagation and induced fluid flow inside the longitudinal fracture was included without sacrificing the computational accuracy and efficiency. Comparing the dispersion and attenuation curves with those from an analytical approach validates the effectiveness of the proposed method. Because the differences in the dispersion and attenuation properties between numerical and analytical results are pretty small, the comparisons clearly showed that the proposed method could accurately simulate 3D elastic wave propagation, 2D fluid flow, and their interaction. The application of the method to the monopole sonic logging situation showed good agreement with the previous study and extended the application range to smaller fracture apertures. We found that the attenuation characteristics are sensitive to small opening fractures and valuable for estimating the fracture aperture and fluid viscosity. Numerical experiments with various fracture apertures and fluid viscosity demonstrated the effectiveness of estimating the aperture and viscosity by analyzing the behavior of the Stoneley waves. Because the method can estimate the aperture and viscosity separately, it could be used for evaluating the efficiency of hydraulic fracturing. The application of the proposed method to noisy data also showed the possibility of the estimation. In the case of a comparatively low S/N (e.g., equal to two), the estimation may include some error because of the disturbance of the attenuation curves.

We thank the associate editor, Y. Liu, and three anonymous reviewers for their thoughtful comments and suggestions, which improved our manuscript significantly.

Data associated with this research are confidential and cannot be released.

APPENDIX A 3D SEISMIC WAVE MODELING

Takekawa et al. (2014a, 2014b) obtain the deformation gradient tensor and the Hamiltonian equation:

Fi=(jrijrij0)Ai1wij,
(A-1)
mjvjt=iFiSiAi1rij0ΔBiwij,
(A-2)
where
Ai=jrij0rij0wij,
(A-3)
and Fi is the deformation gradient tensor, rij is the distance between particles, rij0 is the initial distance between particles, mj is the mass of the velocity particle, Si is the stress tensor, ΔBi is the volume of the stress particle, and wij is the weighting coefficient depending on the distance between particles i and j.
This study also introduces the staggered particles with nonuniform intervals (Figure A-1). Particle spacing (element size) and physical properties are defined at stress particle i. The stress particle is located in the middle of the element, and the displacement (velocity) particle is located at the edges of the element. Each stress particle is surrounded by eight velocity particles. Unit direction vectors toward each velocity particle in the x-, y-, and z-directions are defined as follows:
{lx(18)=(+1,1,+1,1,+1,1,+1,1)ly(18)=(+1,+1,1,1,+1,+1,1,1)lz(18)=(+1,+1,+1,+1,1,1,1,1).
(A-4)
The strain tensor at stress particles can be calculated as follows:
rij0=xi0xj0=(Δxi2lx(j)Δyi2ly(j)Δzi2lz(j)),
(A-5)
Ai=j=18(rij0rij0),wij=j=18wij(Δxi24lx(j)2ΔxiΔyi4lx(j)ly(j)ΔxiΔzi4lx(j)lz(j)ΔxiΔyi4lx(j)ly(j)Δyi24ly(j)2ΔyiΔzi4ly(j)lz(j)ΔxiΔzi4lx(j)lz(j)ΔyiΔzi4ly(j)lz(j)Δzi24lz(j)2)=2wij(Δxi2000Δyi2000Δzi2)(wherewij:const),
(A-6)
Ai1=12wij(1/Δxi20001/Δyi20001/Δzi2),
(A-7)
Fi=j=18(rijrij0)Ai1wij=j=18((rij0+uj)rij0)Ai1wij=I+j=18(ujrij0)Ai1wij=I+j=18uj{(Ai1)Trij0}wij=I+j=18uj(14Δxilx(j)14Δyily(j)14Δzilz(j))=I+fi,
(A-8)
and
Ei=FiTFiI2=(I+fi)T(I+fi)I2=fiTfi+fiT+fi2fiT+fi2,
(A-9)
where Δxi,Δyi, Δzi are the particle spacing in the x-, y-, and z-directions, respectively; Ei is the strain tensor; and fi=FiI, uj is a displacement vector. This expression is consistent with the infinitesimal displacement theory. The stress tensor changing rate is as follows:
Sit=Ci:Eit=Ci:(fitT+fit2),fit=j=18vj(14Δxilx(j)14Δyily(j)14Δzilz(j)),
(A-10)
where Ci is an elastic tensor. Considering that the mass of the velocity particle is calculated by the arithmetic average of the surrounding stress particles, the Hamiltonian equation is given in the following equation:
(i=18ρiΔBi8)vjt=i=18FiSiΔBi(14Δxilx(i)14Δyily(i)14Δzilz(i))i=18SiΔBi(14Δxilx(i)14Δyily(i)14Δzilz(i)),
(A-11)
where ρi is the density. The stress Si and velocity vj are updated at each time step. This algorithm matches the rotated staggered-grid finite-difference method (Saenger et al., 2000) at equally spaced particles. The stability condition is vmaxΔt/Δxmin<1. An adaptive particle arrangement is adopted in this study to minimize the numerical dispersion induced by a coarse particle interval around the borehole. A very fine particle spacing (1 mm interval in the horizontal direction) is used in the vicinity of the borehole. Then, the interval is gradually increased with distance from the borehole. This strategy improves numerical accuracy and computational efficiency (Takekawa et al., 2014c). An absorbing boundary condition (Cerjan et al., 1985) is applied to 2D and 3D wave calculations to suppress artificial reflections from the model edges.

Figure A-2 shows the dispersion property and waveforms obtained by the proposed HPM and analytical methods. The dispersion curve of the analytical solution is found by the grid search algorithm of the dispersion equation, and the analytical waveforms are calculated by the discrete wavenumber integration method (Cheng and Toksöz, 1981; Tang and Cheng, 2004). The comparison result shows good agreement, although discretization brings a small dispersion error. The maximum difference in the slowness between numerical and analytical dispersion curves is approximately 0.9%. This result indicates that the HPM with a nonuniform particle arrangement can accurately simulate the Stoneley wave propagation.

APPENDIX B 2D FLUID FLOW

First, we derive equations 7 and 8. The second term in equation 6 is defined as ζ:
ζ=ikci1ωρfp,
(B-1)
iωζ=kci1ρfp.
(B-2)
By applying the Fourier transformation to equation B-2, equation B-3 can be obtained:
dζdt=kci1ρfp.
(B-3)
Inserting this result into equation 6 leads to
v=kcrL212  μp+ζ.
(B-4)
Next, we explain the discretization of this equation using HPM. In a 2D flow calculation, discretization is performed using an analogy from the difference equation of the HPM with a nonuniform interval shown in Appendix  A. At the stress particle, the pressure is defined, whereas the flow velocity is defined at the velocity particle. Each stress particle is surrounded by four velocity particles. Unit direction vectors toward each velocity particle in the x- and z-directions are described in the following equations:
{lx(14)=(+1,1,+1,1)lz(14)=(+1,+1,1,1).
(B-5)
The spatial derivative of pressure and flow velocity can be expressed as follows:
(i=14ΔBi4)pj=i=14piΔBi(12Δxilx(i)12Δzilz(i)),
(B-6)
1Kpit=·v=j=14(vxj12Δxilx(j)+vzj12Δzilz(j)),
(B-7)
where K is the bulk modulus. The stability condition of the wave equation is the same as in the 3D case, and the stability condition of the diffusion equation is expressed as follows, with kd as the diffusion coefficient:
kdΔt(1(Δxmin)2+(Δzmin)2)12.
(B-8)

Here, 2D and 3D systems have different stability conditions. We need to use a smaller time step of them for stable calculation. The 2D diffusion equation by this algorithm was verified by comparing it with the exact solution p(r,t)=P0/(4πkdt)exp[r2/(4kdt)]. The numerical conditions of the validation model are the same as those of the sonic logging simulation, i.e., L=50  μm, μ=4mPa·s, ρf=1000  kg/m3, αf=1500  m/s, Kf=ρfαf2, and kd=KfL2/(12μ). The comparison result is shown in Figure B-1. The numerical result has excellent agreement with the exact solution. This result indicates that the preceding calculation strategy can reproduce 2D fluid flow with high accuracy.

APPENDIX C THE DISPERSION EQUATION OF THE FRACTURE MODE WHEN FILLED WITH A VISCOUS FLUID

The viscous fluid flow from Tang and Cheng (1989) and the theory of the elastic fracture mode from Ferrazzini and Aki (1987) are combined to derive the theoretical dispersion equation of the fracture mode, including the viscous fluid indicated in Korneev (2008).

Suppose that a 2D fracture extends in the x-direction and opens in the z-direction, with z = 0 at the center of the opening. The displacement potentials ϕ and ψ of the P and S waves satisfy the following wave equations:
2ϕ+ω2VP2ϕ=0,
(C-1)
2ψ+ω2VS2ψ=0,
(C-2)
where VP and VS are the P-wave velocity and S-wave velocity, respectively. Considering only the case where the z-coordinate is positive due to the symmetry, the potential solution can be found as follows:
ϕ=Aepzcos(kx)eiωt,
(C-3)
ψ=Beszsin(kx)eiωt,
(C-4)
where k is the wavenumber in the x-direction and,
p=k2ω2VP2,s=k2ω2VS2.
(C-5)
In this case, the displacement field can be calculated as follows:
ux=ϕxψz=(AB)·(kepzsesz)sin(kx)eiωt,
(C-6)
uz=ϕz+ψx=(AB)·(pepzkesz)cos(kx)eiωt.
(C-7)
In addition, the stress field can be calculated as follows:
τzz=λuxx+(λ+2μ)uzz=(AB)·((λk2+(λ+2μ)p2)epz2μskesz)cos(kx)eiωt,
(C-8)
τzx=μ(uxz+uzx)=(AB)·μ(2pkepz(k2s2)esz)sin(kx)eiωt,
(C-9)
where λ=ρ(VP22VS2) and μ=ρVS2 are the Lamé constants.
Next, in the viscous fluid in the fracture, following Tang and Cheng (1989), we transform the Navier-Stokes equations and the continuity equation into the frequency domain and consider only first-order perturbations. When decomposing displacement into potentials ϕf and ψf for acoustic P waves and viscous shear, the following equations are obtained:
2ϕ+ω2VPf2ϕ=0,
(C-10)
2ψ+ω2VSf2ψ=0.
(C-11)
Here, the complex velocities VPf and VSf of the P and S waves of the fluid are shown as follows:
VPf=αf243iων,VSf=iων,
(C-12)
where αf is the acoustic P-wave velocity and ν is the kinematic viscosity of the fluid.
The potential solution to the fracture can be found using the cosine functions as follows:
ϕf=Afcos(fz)cos(kx)eiωt,
(C-13)
ψf=Bfsin(f¯z)sin(kx)eiωt,
(C-14)
where
f2=ω2VPf2k2,f¯2=ω2VSf2k2.
(C-15)
At this time, the velocity field can be calculated as follows:
uxf=ϕfxψfz=(AfBf)·(kcos(fz)fcos(fz))sin(kx)eiωt,
(C-16)
uzf=ϕfz+ψfx=(AfBf)·(fsin(fz)ksin(fz))cos(kx)eiωt.
(C-17)
In addition, the stress field can be calculated as follows:
τzzf=λfuxfx+(λf+2μf)uzfz=(AfBf)·((λfk2(λf+2μf)f2)cos(fz)2μffkcos(fz))cos(kx)eiωt,
(C-18)
τzxf=μf(uxfz+uzfx)=(AfBf)·μf(2fksin(fz)(k2+f2)sin(fz))sin(kx)eiωt,
(C-19)
where λf=ρf(VPf22VSf2) and μf=ρfVsf2 are the complex Lamé constants and ρf is the density of the fluid. By setting the fracture aperture as L and imposing the boundary condition that the displacement (two components) and stress (two components) are continuous at the formation-fluid boundary (z=L/2), the following equation is obtained:
Gx=0,x=[A,B,Af,Bf]T.
(C-20)
From the conditions for this equation to have a nontrivial solution, we obtain the following dispersion equation (components common to the rows or columns are dropped):
detG=0,
(C-21)
{G11=kG12=sG13=kcos(fL/2)G14=f¯cos(f¯L/2),{G21=pG22=kG23=fsin(fL/2)G24=ksin(f¯L/2),{G31=λk2+(λ+2μ)p2G32=2μskG33=(λfk2(λf+2μf)f2)cos(fL/2)G34=2μff¯kcos(f¯L/2),{G41=2μpkG42=μ(k2s2)G43=2μffksin(fL/2)G44=μf(k2+f¯2)sin(f¯L/2).
(C-22)
The complex wavenumber satisfying the dispersion equation can be estimated, and the analytical solutions of the dispersion and attenuation curves of the fracture mode can be calculated. The complex wavenumber solution is estimated when |detG| takes a minimum value.

The structure of the matrix G is shown in Figure C-1. Imposing only the condition that the wall displacement of the fluid part is zero (determinant = 0 for the submatrix shown in red in Figure C-1), we obtain the characteristic equation 3 for viscous fluid flow in a rigid fracture shown in Tang and Cheng (1989). In addition, by imposing only the condition without viscosity in this matrix (determinant = 0 for the submatrix shown in green in Figure C-1), we obtain the dispersion equation 29 for the elastic fracture mode shown in Ferrazzini and Aki (1987).

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.