## ABSTRACT

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.

## INTRODUCTION

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 10^{12}–10^{13} 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.

## METHODOLOGY

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.

### 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.

### 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 $\Delta x=\Delta z=6.25\u2009\u2009mm,\Delta y=1\u2009\u2009mm$ 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\u2009\u2009\mu m$ is much smaller than the grid spacing $\Delta y$ (1/20). In this case, we can use the analytical dispersion equation of the fracture mode (Ferrazzini and Aki, 1987) as follows:

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:

$\zeta $ is updated in $dt3D$ (not dependent on mixing timing),

$\zeta $ is updated in $dt2D$, instant mixing in $dt2D$,

$\zeta $ is updated in $dt2D$, instant mixing in $dt3D$,

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\xd710\u22129\u2009\u2009s$ for the diffusion calculation in the 2D system and $2.5\xd710\u22127\u2009\u2009s$ 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 $\mu =4\u2009\u2009mPa\xb7s$. Numerical experiments are carried out with the following three time-setting patterns:

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

$dt2D=2.0\xd710\u22129\u2009\u2009s=dt3D/125$, $dp$ is mixed immediately in $dt2D$,

$dt2D=2.0\xd710\u22129\u2009\u2009s=dt3D/125$, $dp$ is mixed in $dt3D$,

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\u2009\u2009\mu m$, $\mu =4\u2009\u2009mPa\xb7s$) 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\xd710\u22129\u2009\u2009s$ for the diffusion calculation and mixing timing of the 2D system and $2.5\xd710\u22127\u2009\u2009s$ for the propagation calculation of the 2D system and seismic wave propagation of the 3D system.

## MODEL

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 $\Delta x=\Delta y=1\u2009\u2009mm$, $\Delta z=6.25\u2009\u2009mm$ 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.

## RESULT

Figure 11 shows an example of the simulation result in the propagation-dominant condition ($L=50\u2009\u2009\mu 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\u2009\u2009\mu m$ begins to move away from that without fracture at approximately 1 kHz, whereas the slowness with the aperture of $100\u2009\u2009\mu 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\u2009\u2009\mu m$, $\mu =4\u2009\u2009mPa\xb7s$. 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\u2009\u2009\mu m$, $\mu =4\u2009\u2009mPa\xb7s$. 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\u2009\u2009\mu m$, $\mu =4\u2009\u2009mPa\xb7s$, 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.

## ESTIMATION OF APERTURE AND 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.

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\u2009\u2009\mu 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.

## CONCLUSION

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.

## ACKNOWLEDGMENT

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

## DATA AND MATERIALS AVAILABILITY

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:

*x*-,

*y*-, and

*z*-directions are defined as follows:

*x*-,

*y*-, and

*z*-directions, respectively; $Ei$ is the strain tensor; and $fi=Fi\u2212I$, $uj$ is a displacement vector. This expression is consistent with the infinitesimal displacement theory. The stress tensor changing rate is as follows:

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

*x*- and

*z*-directions are described in the following equations:

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\pi kdt)exp[\u2212r2/(4kdt)]$. The numerical conditions of the validation model are the same as those of the sonic logging simulation, i.e., $L=50\u2009\u2009\mu m$, $\mu =4\u2009mPa\xb7s$, $\rho f=1000\u2009\u2009kg/m3$, $\alpha f=1500\u2009\u2009m/s$, $Kf=\rho f\alpha f2$, and $kd=KfL2/(12\mu )$. 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).

*x*-direction and opens in the

*z*-direction, with

*z*= 0 at the center of the opening. The displacement potentials $\varphi $ and $\psi $ of the P and S waves satisfy the following wave equations:

*z*-coordinate is positive due to the symmetry, the potential solution can be found as follows:

*x*-direction and,

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.