Based on the mechanism of the interaction between elastic wave and fluid transport property of natural or artificial fractures, the microseism, vertical seismic profile (VSP) and tube wave measurement have already been used in the fracture characterization ^[1,2,3]. However, as the oil and gas exploration is extending to the “deep sea and deep formation” in China, the increasing detection depth raises the cost and reduces the measurement accuracy, consequently influencing the evaluation effectiveness of these methods.

Acoustic logging is a common downhole detection method, with advantages of in-situ measurement and high accuracy. The dynamic fluid flow induced by the low frequency Stoneley-wave during acoustic logging in a permeable formation has been well studied theoretically ^[4-5] and experimentally ^[6,7], yielding many practical methods to characterize the fractures using Stoneley-wave reflection, transmission, attenuation and dispersion ^[8,9]. However, these methods, which are mainly used for natural fracture evaluation, have rarely been used in hydraulic fracture evaluation. Hydraulic fractures are known to have fractal distribution characteristics ^[10], and their fluid transport property varies dramatically. The fracture width is of hundred-micron level (the permeability is hundreds of square-microns) when the fracture is open and well-propped, but reduces to the micron-level or even less (the permeability is several to dozens of square-microns) after fracture closure. The change of fracture width leads to a complicated change in the acoustic field around the fracture. Another fact is that the cementation and perforation of the borehole are commonly required before hydraulic fracturing, so the fracture fluid flow due to the borehole acoustic field can only induced through the perforation holes. Therefore, the fluid coupling in a cased borehole should be significantly different from that in an open borehole and needs to be modified.

The specific factors discussed above are fully studied in this study. We first built the permeability model of a hydraulically fractured formation based on the laboratory experiment samples, and then solved the diffusion-type wave equation of dynamic fluid flow using the finite-difference method. The dynamic fluid flow in hydraulic fractures is then coupled with the borehole acoustic field through perforation holes using the simplified Biot-Rosenbaum (B-R) method. We thoroughly analyzed the relationship between Stoneley-wave velocity dispersion behaviors and fracture permeability, proposed an evaluation method for the fluid transport property of hydraulic fractures. Finally, this novel method has been validated with field data analysis.

The hydraulically fractured formation is a complex and inhomogeneous permeable medium consisted of the low permeable bedrock and hydraulic fractures. With a given angular frequency ω, the dynamic fluid flow at x inside the medium is described by the diffusion-type wave equation ^[11]:

In Eq. (1), α is the dynamic fluid diffusivity, which is described by:

In Eq. (2), ξ is the correction factor of the medium framework elasticity and equals to 0 assuming that bedrock is incompressible compared to fluid. K (ω, x) is the dynamic permeability given by Johnson et al. ^[12]:

In Eq. (3), K₀(x) is the static permeability at ω=0, which varies with spatial coordinate at x;τ and m are relevant to the shape of pores or fractures ^[12]: for τ≈3 and m=2, ϕ represents the actual porosity in a typical porous medium; while for τ=1 and m=3, ϕ is taken as 1 to represent a fracture-infill medium. The variation of K (ω, x) with the frequency dictates that the fluid flow field governed by Eq. (1) is the viscous flow at low frequencies, and becomes the propagational wave at high frequencies ^[11].

K₀(x) in Eq. (3) varies spatially because of the heterogeneity of hydraulically fractured formation. We assume that hydraulic fractures are similar along the borehole axis (z axis), only the variation on the (xOy) plane perpendicular to the borehole axis needs to be analyzed, that is, K₀(x)=K₀(x, y), simplifying the fractured formation simulation to a 2-dimensional problem. Many scholars have studied the hydraulic fracture distribution through experiments. These experiments usually used a reduced scale cement (or rock) as base rock with a perforated steel pipe to simulate a perforated borehole. Loading the specimen with confining pressure, the high-pressure fluid was injected into the perforated pipe to induce fractures. Lee et al. ^[10] have shown that these stimulated fractures network ^[13] have fractal characteristics, and can closely reflect the shape of hydraulic fractures. For the model shown in Fig. 1a, we used the threshold segmentation method to acquire the fracture distribution from the specimen ^[13].

The irregular fracture surface results in the permeability variation with spatial position. This study uses a fractal fracture surface modeling method proposed by Brown ^[14]. Taking the minimum (full contact, equals to 0) and maximum fracture width as reference values, the fracture widths are generated according to the fractal distribution inside the hydraulic fracture profile. The overall statistic of generated fracture widths shows a normal distribution, and its average is about half of the maximum. The fracture width d at each position is transferred into the static permeability as K₀=d²/12 ^[15]. Taking the fracture distribution in Fig. 1a as an example, the model is first scaled to the actual borehole dimension, and then the permeability map is generated for the scaled model (Fig. 1b), where the model size and grid size are 4.096 m × 4.096 m and 256×256, respectively. The normalized permeability map is shown in Fig. 1b, where the values of 1 and 0 indicate the maximum and minimum permeability in this model, respectively, with the variation in between showing the complicated spatial distribution.

The fractal permeability distribution model based on the hydraulic fracturing experiment lays the foundation for numerical simulation. The bedrock (tight sandstone) permeability and porosity are set as 1×10^-3 μm² and 8% respectively ^[16]; the elastic parameters of bedrock ^[17] and cased borehole are given in Table 1. The same viscosity value of 1.14 mPa·s is used for both borehole and formation fluid.

During hydraulic fracturing, the injection fluid carries proppants (e.g. quartz sand) into the formation. In order to maximize the stimulated volume and optimize the propping effect, the commonly used proppant diameter is usually 0.15-2.00 mm (100-10 mesh). Under the action of formation stress, the fracture has undergone multiple stages of initial propping and opening, mid-term compaction, and final steady-state. The fracture width changes from the initial hundreds of microns to the final value of a few microns, and the corresponding permeability also varies from hundreds to tens or even a few square microns.

Eq. (1) is solved for the abovementioned model and parameters. Because of the inhomogeneous variations of permeability distribution, Eq. (1) is solved using the ADI (Alternating Direction Implicit) finite difference method because of its good numerical stability ^[18].

Excited by the Stoneley-wave propagation along cased borehole with the wavenumber k_e, the formation fluid pressure in Eq. (1) has the form can be written as $p{{\text{e}}^{i{{k}_{\text{e}}}z}}$, where z is the borehole distance, Eq. (1) can then be rewritten as:

For Eq. (4), the boundary conditions inside and far away from the borehole are:

With the boundary conditions, the dynamic fluid pressure (fluid pressure for short) p(x, y) at a given frequency can be obtained by solving Eq. (4) using ADI finite difference method (see solution details in reference [18]).

For hydraulically fractures with the maximum permeability (K_max) of 200 μm² (equivalent to the maximum fracture width of 50 μm), the average fracture width and corresponding average permeability (${K}'$) are 25 μm and 50 μm², respectively. Because of the equivalency of the fracture width and permeability, only the permeability will be discussed in the later text. For this fracture model, the fluid pressure field at the frequency of 500 Hz is shown in Fig. 2a, where the pressure amplitude is nor-malized by its maximum value and displayed using a color scheme, with brighter color representing higher fluid pressure and vice versa. In comparison, the fluid pressure for the plane fracture model with the permeability of 50 μm² is given in Fig. 2b. The comparison shows that the dynamic fluid flow in the hydraulic fracture is confined inside the fracture; whereas, the contribution of bedrock to dynamic fluid flow is negligible due to the low background permeability.

The fluid pressure along a profile in both models of Fig. 2 (see blue dashed lines) is shown in Fig. 3. The colors of black, red and blue indicate the amplitude, real part and imaginary part of fluid pressure, respectively. The fluid pressure field of hydraulically fractured formation concentrates in the range of 0-2 m around the borehole, which marks the sensitivity range of the 500 Hz borehole acoustic wave. The dynamic fluid flow is locally controlled by the complex distribution and irregular shape of fracture surface, thus cannot be directly compared with that of a plane fracture model. However, it is the overall effect of fluid diffusion in the fracture that interacts with the borehole acoustic wave, and this effect can be examined by the analysis of the borehole acoustic field.

The B-R theory ^[19] describes the interaction between borehole acoustic field and a permeable formation. The simplified B-R theory is an analytical expression for the low-frequency Stoneley-wave application of the theory ^[11], and has been validated theoretically and experimentally. The simplified B-R theory ^[11] has been effectively applied to field acoustic data to obtain the formation permeability. However, the application of the theory to the hydraulically fractured cased borehole needs to consider the effects of perforation holes, as shown in Fig. 4, where the perforated tunnels inside formation are neglected assuming they are destroyed after hydraulic fracturing.

Including the effect of casing perforation in the simplified B-R theory and following the theoretical derivation from Tang et al. ^[11], the borehole Stoneley-wave wavenumber is modified from the original theory to become:

Essentially, the theory combines the elastic wavenumber k_e and the influence of dynamic fluid flow on the borehole wall, which is the second term under the square root in Eq. (6), where β is the perforation factor relevant to the perforated hole area, perforation density and borehole radius etc. β ranges between 0 to 0.2 according to common perforation parameters and is set as 0.1 in this study; R is the inner radius of casing, which equals to the borehole fluid outer radius in Table 1; a is the outer radius of the acoustic logging tool, and usually less than 0.05 m; ${{\left. \frac{{{U}_{\text{f}}}}{p} \right|}_{r=R}}$, known as the borehole wall conductance, is the ratio of dynamic fluid flux and the fluid pressure at the borehole wall. For a homogeneous formation, the wall conductance can be calculated using Eq. (7) ^[4]:

The wall conductance of a hydraulically fractured formation is much more complicated, as it must be numerically calculated using the finite-difference result of dynamic fluid flow. For this calculation, the result of Eq. (4), p(x, y) calculated in the Cartesian coordinates, must be converted into the polar coordinates using Eq. (8). In the conversion, the fluid pressure values on the Cartesian grids are interpolated onto the polar coordinates.

Then the average conductance at the wall of full borehole is obtained by substituting p(r,θ) into Eq. (9) ^[18]:

The radial derivative and azimuthal integration in Eq. (9) are numerically calculated by difference and summation, respectively. The wall conductance is substituted into Eq. (6) to get the Stoneley-wave wavenumber, and then the Stoneley-wave phase velocity dispersion (Stoneley dispersion for short) is calculated using Eq. (10):

In Eq. (10), Re denotes taking the real part of variable (the same in the following).

The wall conductance of the dynamic fluid flow directly reflects the fluid transport property of the hydraulically fractured formation, which can be used to validate the numerical results.

Taking models (${K}'\text{=50 }\mu {{\text{m}}^{\text{2}}}$) in Fig. 2 as an example, the wall conductance and Stoneley dispersion values of hydraulic fracture and plane fracture models calculated by analytical solutions and numerical solutions in the frequency range of 0 to 5 kHz are shown in Fig. 5a and 5b, respectively. The wall conductance values of the plane fracture model obtained by Eqs. (7) and (9) are consistent (Fig. 5a), and the analytical and numerical solutions almost coincide with each other. The Stoneley dispersion calculated by substituting them into Eqs. (6) and (10) is also in good agreement (Fig. 5b), which further verifies the accuracy of the numerical solution. However, the wall conductance and Stoneley dispersion of the hydraulic fracture model both deviate from their plane fracture model counterpart, particularly in the low frequency range of 0 to 2 kHz. On the one hand, this difference comes from the strong inhomogeneous changes of fracture morphology and fracture width; on the other hand, the quadratic relationship between permeability and fracture width causes the average permeability of fractures to be slightly larger than the average fracture width permeability. Using the average permeability ($\bar{K}\approx 70\text{ }\mu {{\text{m}}^{\text{2}}}$) instead of the average fracture permeability (${K}'\text{=50 }\mu {{\text{m}}^{\text{2}}}$), the new analytical result using Eq. (7) is more in line with the numerical result of the hydraulically fractured model (Stoneley dispersion curves of hydraulic fracture and equivalent plane fracture in Fig. 5b). Thus, it can be considered that the Stoneley dispersion curve reflects the average permeability of the hydraulic fractures in a hydraulically fractured formation. To further investigate the Stoneley dispersion variation with the average permeability change due to different fracture widths, the Stoneley dispersion curves with a variety of average permeability values (25, 50, 75, 100 μm²) are calculated using Eqs. (6) and (7), as shown in Fig. 6. The other parameters in the calculations are listed in Table 1, and the tool effect is neglected (a=0). Because the B-R theory is commonly used in low and medium permeability formations, it is important to discuss its application results for highly permeable hydraulically fractured formations.

For the low permeability (25 μm²) scenario of Fig. 6, Stoneley wave velocity decreases monotonically with the decrease of frequency. This decrease trend becomes more prominent with the increase of permeability (50 μm²). Interestingly, as the permeability further increases to 75 μm² and 100 μm², the Stoneley-wave velocity increases to a peak value in a certain frequency band. The wave velocity drops rapidly below this frequency band, but decreases slowly to approach a constant above this frequency band. This phenomenon can be explained by the diffusion-type wave equation (Equation (1)). Under dynamic stimulation, dynamic fluid flow is limited to a viscous skin depth range on the fracture interface. For example, when the frequency is 1 kHz, the viscous skin depth of fracture saturated with water is 20 μm. The fluid diffuses within the viscous skin depth. Outside the skin depth, the fluid motion is characterized as wave propagation. At low frequency, the viscosity skin depth is large, and the fluid motion governed by Equation (1) is mainly viscous diffustion, and its coupling with the borehole acoustic field reduces the stoneley wave velocity. However, when fracture permeability and frequency increase (viscous skin depth is negligible), Eq. (1) is transformed to a wave equation, and the dynamic fluid flow in the fracture is analogues to the propagational wave in a plane fracture, also known as the Krauklis-wave, as studied by Krauklis ^[20], Ferrazzini and Aki ^[21]. The absence of viscous effect at the borehole wall reduces the coupling of the borehole with the acoustic field, such that the Stoneley-wave velocity tends to approach the free-space acoustic velocity. As frequency further increases, the Stoneley velocity drops slightly, and finally tends to the Scholte wave velocity, which is a little lower than the fluid velocity ^[22].

The average permeability of hydraulic fractures and the analytical result of Eq. (7) provide a fast and effective forward calculation for the following evaluation method. Meanwhile, the theoretical and numerical analyses indicate an important phenomenon that the high permeability (related to the open hydraulic fracture cases) produces a Stoneley dispersion peak around 1 kHz. Such phenomenon has not been reported in conventional permeable reservoir ((0.1-2000) × 10^-3 μm²) logging, and is a special characteristic for hydraulically fractured formation ^[16]. It can be regarded as an indicator of the hydraulic fracturing effectiveness.

The previous simulation analyzed the Stoneley dispersion characteristics of the equivalent average permeability model for hydraulic fractures extending along the borehole axis. In many cases, e.g., horizontal wells with dominant vertical stress, the opening height of the hydraulic fractures at each perforation hole position may be smaller than the perforation spacing, i.e., the fracture cluster is not well-connected along the borehole axis direction. We therefore propose the general fracture cluster model in Fig. 7a, in which the fractured layers extending from perforation holes are alternatively distributed along the borehole axis. Using the equivalent plane fracture to represent the hydraulic fracture (Fig. 5), the model can be further simplified as fracture cluster model in Fig. 7b. Between each two fractured layers in this model, the fractured layer thickness is h_frac and the unfractured formation thickness is h_fm. This defines the hydraulic fracture connectivity along the borehole axis, γ=h_frac/(h_frac+ h_fm), for assessing the effectiveness of hydraulic fracuring. Because (h_frac+h_fm) is the perforation spacing (a parameter related to perforation density), it is obvious that γ[0,1]. Particularly, if the fractured layers are well connected, i.e. the hydraulic fractures are fully opened along the borehole axis (h_fm=0), one then has γ=1; on the other hand, if the formation is not fractured (h_frac=0), one has γ=0. Note that the perforation spacing usually ranges from 0.02 m to 0.06 m in most perforation guns, h_frac+h_fm is therefore far less than the Stoneley-wave wavelength (greater than 1 m at 1 kHz). Under this condition, the equivalent medium averaging method can be used to get the Stoneley- wave velocity as the weighted average over the fractured and unfractured formation intervals:

where the weight is γ (0<γ<1); k_fm is the unfractured formation Stoneley wavenumber, i.e., the k_e parameter in Eq. (6); the fractured formation Stoneley-wave wavenumber k_frac is calculated using Eqs. (6) and (7). Eq. (11) provides an effective evaluation model to characterize the fluid transport property of stimulated hydraulic fractures for both vertical and horizontal wells.

Fig. 8 shows the Stoneley-wave velocity dispersion of Eq. (11) for two average fracture permeabilities of 25 (Fig. 8a) and 100 (Fig. 8b) μm² and for various connectivity values from 0 to 1, the dispersion curve shows the highest velocity reduction (increase) when the fractures are connected (γ=1). With decreasing connectivity (γ<1), the effect of fracture deminishes. When the fracture connectivity goes to 0, the effect of fracture disappears and the Stoneley dispersion curve approaches its elastic-formation counterpart. The modeling examples show that the fracture connectivity is an important evaluation parameter, which reflects the connection state of hydraulic fracture clusters and the effectiveness of the fracturing operation.

The field data analysis workflow using our evaluation method is described as follows. We first extract the Stoneley dispersion data from the acoustic logging of a hydraulically fractured formation. We then obtain other parameters, such as fluid and tool parameters, perforation factor and calculate the elastic Stoneley-wave wavenumber. These parameters are either obtained from the logging data, or calibrated from a non-fractured interval. The dis- persion extraction and parameter calibration methods are described in reference [22]; the perforation factor is obtained from the perforating tool. Finally, the following objective function is constructed to obtain the average permeability and the connectivity of the hydraulic fractures through an inversion procedure.

In Eq. (12), v_m is the theoretical Stoneley-wave velocity calculated using Eq. (6), and f is frequency. In the frequency range of f_min to f_max, $\bar{K}$(average permeability) and γ (hydraulic fracture connectivity) of the hydraulic fractures are obtained by minimizing the objective function F.

The example from a horizontal production well is given in Fig. 9. The hydraulic fracturing was conducted at two intervals of 1183-1187 m and 1192-1194 m, respectively. The gamma-ray and P-wave slowness curves before and after fracturing show very little difference (Fig. 9a), indicating no overall changes in the formation. The P-wave radial travel time tomography technique ^[23] is used to obtain the P-wave radial velocity profile after fracturing (Fig. 9a). This technique can assess the radial velocity change of the formation around the borehole. The color variation indicates the velocity changes at the fractured interval, whereas in the unfractured formation the radial velocity change is minimal. The mechanism is that the fracturing stimulated fractures around the borehole reduce the wave velocity. More importantly, it is noted that the Stoneley-wave (center frequency around 1.5 kHz) has an obvious “onward trend” in the fractured interval comparing to the unfractured one, which means the Stoneley-wave in the fractured interval experiences a velocity increase. The velocity increase is closely related to the low-frequency Stoneley-wave velocity peak for a high-permeability fractured formation.

By comparing the Stoneley-wave array waveforms at the midpoint (1185 m) of the upper fractured interval and the reference point (1177 m) of the unfractured interval (Fig. 9b), we found an onward trend of waveforms in the fractured interval comparing to the unfractured one, as shown by the two crests with the arrows pointed. The extracted Stoneley dispersion data also show significant differences. In Fig. 9c, the Stoneley-wave velocity of the unfractured interval slightly varies with the frequency, whereas there is an obvious maximum at 1 kHz in the fractured interval case. The dispersion data are inverted using Eq. (12). For the inversion, the velocity and viscosity of fracture fluid are given at 1330 m/s and 1.1 mPa·s, respectively, and the other parameters are obtained from logging data or roughly estimated. The inversion gives the hydraulic fracture connectivity value of 0.48 and the average permeability value of 330 μm². This order of permeability, according to conductivity test results of API (American Petroleum Institute), corresponds to the permeability value of well-propped open fractures ^[24]. The inversion-fitted Stoneley dispersion curve (Fig. 9c) is in a good agreement with the measured dispersion data of the fractured interval (red), both showing a prominent velocity peak around 1 kHz, validating the theoretical prediction for the highly permeable fracture case. This example also further confirms that the onward trend of the waveform and the velocity peak in the dispersion data in a fractured interval can be taken as an effective indicator of well-developed and highly permeable hydraulic fractures.

Fracturing evaluation also needs to determine the changes of hydraulic fractures. During the production, the fluid transport property of hydraulic fractures may be reduced because of proppant deformation, formation pressure change, and blockage etc. For the vertical well example in Fig. 10, the hydraulic fracturing was conducted at the interval of 3963-3969 m. In Fig. 10a, the attenuation (brown envelope) and anisotropy value of S-wave before and after fracturing has confirmed the hydraulic fracture development around the borehole ^[10]. However, the Stoneley-wave (with a center frequency about 1.5 kHz) shows a delay in the fractured interval, which implies a velocity drop, instead of increase, as compared to the waveform “onward” and dispersion curve ‘’peaking’’ phenomena in Fig. 9.

Similarly, we compared the Stoneley-wave array waveforms recorded in the fractured interval (3966 m) to their counterpart recorded in the unfractured interval (3979 m) (Fig. 10b). The comparison shows a clear waveform delay for the fractured interval waveform, as indicated by the relative shift beween the two sets of waveforms. The extracted Stoneley dispersion data also show significant differences. In Fig. 10c, the overall Stoneley dispersion of the fractured interval is lower than that of the unfractured interval, indicating an overall velocity reduction over the frequency range. The inversion for this case gives the average permeability of 15 μm² and the fracture connectivity of 1, respectively by the inversion-fitted theoretical Stoneley dispersion curve fits well with the data. There are two possible explanations for this phenomenon. The first one is that the hydraulic fractures are not fully developed after fracturing. The other one is that these fractures are undergoing gradual closure after fracturing, and, with the fracture closure, fracture permeability decreases. Looking from the significant S-wave anisotropy (caused by hydraulic fractures) after fracturing, the latter explanation is more reasonable.

The Stoneley-wave velocity dispersion characteristics caused by complex hydraulic fractures can be equivalent to those of the plane fracture model, provided that the average permeability of the hydraulic fracture model is equated to the permeability of the plane fracture. Open and well-developed hydraulic fractures can increase the Stoneley-wave velocity and lead to the “onward trending” phenomenon of Stoneley-wave (especially the arrival portion) at low frequencies. In contrast, the ineffective fracturing operation or the fracture closure during the production results in the overall Stoneley-wave velocity reduction, producing a “delay trending” phenomenon on the waveform. The effects of both scenarios have been observed in field data analyses and can be used to indicate fracturing effectiveness.

Based on the theoretical modeling and field data analysis results, a fluid transport property evaluation method for hydraulic fractures has been proposed to estimate the average permeability and hydraulic fracture connectivity of the fractures. The field data analysis results have proved the feasibility and validity of this evaluation method.

a—outer radius of the tool, m;

d—fracture width, m;

f—frequency, Hz;

f_min, f_max—upper and lower frequency limits of velocity dispersion extraction, Hz;

F—fitting objective function, m²/s²;

h_fm—unfractured formation layer thickness, m;

h_frac—fractured layer thickness, m;

k—Stoneley-wave wavenumber, m^-1;

k_e—elastic Stoneley-wave wavenumber, m^-1;

k_fm—Stoneley-wave wavenumber of the unfractured formation layer, m^-1;

k_frac—Stoneley-wave wavenumber of the fractured layer, m^-1;

K—dynamic permeability, m²;

K_max—maximum (static) permeability of the hydraulic fracture, m²;

${K}'$—(static) permeability corresponding to the average fracture width, μm²;

$\overline{K}$—average (static) permeability, μm²;

K₀—static permeability, m²;

K_f—fluid modulus, Pa;

K₀, K₁—0^th and 1^st of second kind of modified Bessel function;

p—fluid pressure, Pa;

p₀—Stoneley-wave amplitude (sound pressure), Pa;

r—radial coordinate in the polar coordinates, m;

R—inner radius of the casing, m;

U_f—fluid flux at the borehole wall, m;

${{\bar{U}}_{\text{f}}}$—average dynamic fluid flux at the borehole wall, m;

v—equivalent Stoneley-wave velocity, m/s;

v_m, v_d—Stoneley-wave velocity dispersion extracted from theoretical calculation and the field data, m/s;

v_st—Stoneley-wave phase velocity, m/s;

x, y, z—Cartesian coordinates, m;

x—spatial vector, m;

α—dynamic fluid diffusivity, m²/s;

β—perforation factor;

γ—hydraulic fracture connectivity;

θ—angular coordinate in the polar coordinates, rad;

μ—fluid viscosity, Pa·s;

ξ—correction factor of medium framework elasticity, dimensionless;

ρ—fluid density, kg/m³;

τ, m—coefficients related to the shape of pore or fracture, dimensionless;

ϕ—porosity, %;

ω—angular frequency, rad/s.