METHOD FOR ESTIMATING PERMEABILITY OF FRACTURED ROCK FORMATIONS FROM INDUCED SLOW FLUID PRESSURE WAVES

Abstract

An embodiment in accordance with the present invention includes a method for estimating the permeability of fractured rock formations from the analysis of a slow fluid pressure wave, which is generated by pressurization of a borehole. Wave propagation in the rock is recorded with TFI™. Poroelastic theory is used to estimate the permeability from the measured wave speed. The present invention offers the opportunity of measuring the reservoir-scale permeability of fractured rock, because the method relies on imaging a wave, which propagates through a large rock volume, on the order of kilometers in size. Traditional methods yield permeability for much smaller rock volumes: well logging tools only measure permeability in the vicinity of a borehole. Pressure transient testing accesses larger rock volumes; however, these volumes are much smaller than for the proposed method, particularly in low-permeability rock formations.

Description

FIELD OF THE INVENTION

The present invention relates generally to geophysical exploration. More particularly, the present invention relates to a method for estimating the permeability of the transmissive fracture network of fractured rock formations from the analysis of a slow fluid pressure wave, which is generated by pressurization of a borehole and imaged with TFI™.

BACKGROUND OF THE INVENTION

Rock permeability is one, if not the most, important parameter to know when it comes to understanding natural and engineered flows in the permeable subsurface. Likewise it is crucial to quantify changes in rock permeability, e.g., due to reservoir stimulation through hydraulic or propellant fracturing. Permeability fields in fractured rock are highly heterogeneous with most of the flow concentrated on transmissive fracture networks. Nonetheless fractures can be rather space-filling, and therefore the transmissive fracture network which predominantly participates in subsurface flow can be rather space filling too. FIG. 1A is a map view of a Tomographic Fracture Image (TFI™) of a space-filling fracture network at the reservoir scale. However, it also possible that there are localized conduits for fluid flow, which can, for example, be linked to fault zones. FIG. 1B illustrates a graphical view of localized zones with relatively high permeability.

Reservoir-scale permeability is an extremely difficult parameter to measure in the deep subsurface. Rock core permeabilities are not representative of the permeabilities of the reservoir-scale fracture networks, which typically exceed the matrix permeability of the cores by orders of magnitude. Pumping tests [e.g., pressure transient analysis (PTA)], in which fluid is injected or withdrawn from a well while changes in borehole pressure are monitored, are not reliable for the low permeabilities often encountered at depth or yield a permeability that is representative only for a small rock volume in the vicinity of the borehole.

Recent advances in rock imaging technology now allow not only mapping the fracture network in the deep subsurface but also identifying dominant flow paths associated with transmissive fracture/fault systems. This is possible through the usage of TFI™, which processes passive seismic data from buried and surface based seismic networks via Seismic Emission Tomography (SET). TFI™ can be used to image directly the dynamic behavior of fracture/fault systems at scales ranging from 10 meters to kilometers. Therefore TFI™ generates five dimensional (5D) semblance data, i.e., the 3D position, time, and seismic energy. TFI™ has, for example, been used to monitor fracture growth during hydraulic fracturing.

Using TFI™, it has been possible to observe slow fluid pressure waves induced by pressurization of a borehole due to the injection of slurries. These so-called Pf waves propagate with a speed on the order of 10 m/s through the transmissive fracture/fault network. The rock domain, which the Pf wave excites, is associated with high-permeability zones within the transmissive fracture network. Due to the small velocity of the Pf waves, it is clear that they are not regular compressional or shear waves, which in natural subsurface environments have velocities on the order of one thousand m/s (e.g., 500 m/s; 1,000 m/s; 2,000 m/s; etc.). So far, a theoretical explanation for the Pf wave has not been provided. Therefore, it is currently not possible to determine rock parameters from the measured speed and attenuation of the Pf wave.

Accordingly, there is a need in the art for a method to estimate the permeability of fractured rock formations from the analysis of a slow fluid pressure wave.

SUMMARY OF THE INVENTION

The foregoing needs are met, to a great extent, by the present invention which provides a method of estimating permeability of fractured rock formations. The method includes generating a Pf wave in a perforated portion of a borehole within the rock formation for the purpose of permeability determination for the fractured rock formation. The method also includes imaging of the Pf wave with TFI™ and analyzing TFI™ data in order to determine the transmissive fracture network, in which the Pf wave propagates as well as the front position of the Pf wave as a function of time within this network. The method includes analyzing the front position data in order to determine the speed of the Pf wave as a function of 3D position over time within the rock formation. The method also includes using poroelastic theory to estimate the 3D permeability field from the 3D field for the Pf wave speed.

In one aspect, the present disclosure provides a method of estimating permeability of a rock formation that includes the steps of generating, by a signal generator, a pressure (Pf) wave in a perforated portion of a borehole within the rock formation; imaging, by an imaging device, the Pf wave to produce imaging data; analyzing, by a controller, the imaging data to identify a front position of the Pf wave as a function of time and to formulate a transmissive fracture network in which the Pf wave propagates; determining the speed of the Pf wave as a function of three dimensional (3D) position within the rock formation based on the front position data; and establishing an estimated 3D permeability field from the 3D field for the Pf wave speed using poroelastic theory.

In one embodiment, the signal generator is borehole pressurization. In one embodiment, the borehole pressurization results from injection of a liquid, a gas, or a supercritical fluid into the borehole. In one embodiment, the liquid is selected from the group consisting of water, isopropyl alcohol, two-butoxyethanol, ethylene glycol, and a slurry.

In one embodiment, the imaging device is a Tomographic Fracture Imaging (“TFI™”) device or a Seismic Emission Tomography (SET) device. In one embodiment, the TFI™ device comprises a surface array or a sub-surface array. In one embodiment, the surface array is a beam forming surface array. In one embodiment, the sub-surface array is positioned 1-10 meters below the surface. In one embodiment, the poroelastic theory is Biot's theory.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings provide visual representations, which will be used to more fully describe the representative embodiments disclosed herein and can be used by those skilled in the art to better understand them and their inherent advantages. In these drawings, like reference numerals identify corresponding elements and:

FIGS. 1A-1B illustrate a map view image of a space-filling fracture network that has been imaged with TFI™ and a graphical view of localized zones with relatively high permeability that could appear in fault zones, respectively.

FIG. 2 shows a series of map view TFI™ images of a Pf wave in the Barnett Shale. The five map views show the semblance for five different times. The gray-scale stands for the magnitude of semblance, with black and white indicating high and low values. The images show that the Pf wave travels about 5,000 ft in 3 minutes, which corresponds roughly to a velocity of 6.3 m/s.

FIGS. 3A-3E illustrate graphical views of how TFI™ data is processed in order to determine the permeability of the transmissive fracture network. FIG. 3A illustrates that first the entire ambient fracture network is mapped by processing passive seismic data with TFI™. Shown here are the central surfaces of the fracture network. FIG. 3B illustrates that the transmissive fracture network is then illuminated by increasing the borehole pressure, which induces a Pf wave. The recorded 5D semblance data is used to determine the arrival times of the Pf wave front in the fractures. Only fractures that form the transmissive fracture network are shown. The color of a fracture indicates the travel time of the Pf wave. FIG. 3C illustrates that from the travel time field, the velocity of the Pf wave is calculated. FIG. 3D illustrates that 5D semblance data from the passive seismic imaging is used to identify the permeable damage zones (shown in gray) around the central surfaces of the fracture network, because these zones also carry fluid flow. This rock volume contains the dominant flow paths and forms the transmissive fracture network. FIG. 3E illustrates using Biot theory, permeability in the transmissive fracture network is estimated from the Pf wave velocity.

FIG. 4 illustrates graphical views of the bottom borehole pressure p_b and flow rate as a function of time as observed during hydraulic fracturing.

FIG. 5 illustrates a graphical view of the dependence of the group velocity on permeability, which was calculated from Biot's theory for poro-elastic media.

DETAILED DESCRIPTION

The presently disclosed subject matter now will be described more fully hereinafter with reference to the accompanying Drawings, in which some, but not all embodiments of the inventions are shown. Like numbers refer to like elements throughout. The presently disclosed subject matter may be embodied in many different forms and should not be construed as limited to the embodiments set forth herein; rather, these embodiments are provided so that this disclosure will satisfy applicable legal requirements. Indeed, many modifications and other embodiments of the presently disclosed subject matter set forth herein will come to mind to one skilled in the art to which the presently disclosed subject matter pertains having the benefit of the teachings presented in the foregoing descriptions and the associated Drawings. Therefore, it is to be understood that the presently disclosed subject matter is not to be limited to the specific embodiments disclosed and that modifications and other embodiments are intended to be included within the scope of the appended claims.

The present invention is based, at least in part, on the discovery that the permeability of reservoir scale (e.g., about 1-10 km in size) geological features and/or transmissive fractures/faults may be estimated by using a signal generator to produce a slow fluid wave that emanates omni-directionally from the signal generator, imaging the slow fluid wave in three dimensions using a geological imaging device/technique (e.g., TFI™), identifying the wave front of the slow fluid wave in three dimensions, calculating the velocity of the wave front in all three dimensions as the slow fluid wave moves outwardly through the geological features and/or transmissive fractures/faults, and analyzing the wave front velocity in three dimensions with a poroelastic theory (e.g., Biot's theory) to calculate permeabilities throughout the geological feature.

Tomographic Fracture Imaging (“TFI™”) is an imaging method that combines Seismic Emission Tomography (SET) with empirical data on fracture geometry to directly image geological fracture networks (e.g. naturally occurring fracture networks, induced fracture networks, and the like). SET is a technique that may identify sources of seismic energy contained within the volume of the geological feature (e.g. a reservoir) being imaged as described in, for example, U.S. Pat. No. 9,389,326 (incorporated herein by reference in its entirety). The seismic energy may be recorded by a sensor array (e.g., a beam forming surface array) and collected and processed to provide a three dimensional (3D) grid of voxels. The semblance value for the seismic energy associated with each voxel may then be calculated and assigned to a node point and then color contoured for semblance value. Advantageously, the seismic data for TFI™ may be collected by a surface array, which enables it to be used to map fracture networks at both the well and reservoir level scale with a high level of resolution. It is also contemplated within the scope of the invention that the seismic data for TFI™ may be collected by a shallow sub-surface sensor array (e.g., a sensor array that has been buried below the surface). In an exemplary embodiment, such a sub-surface sensor array may be positioned 1-100 meters below the surface, or in some cases even deeper.

An embodiment in accordance with the present invention includes a method for estimating the permeability of fractured rock formations from the analysis of a slow fluid pressure wave, which may be generated by pressurization of a borehole. In this regard, borehole pressurization may be used as a signal generator to create the slow fluid pressure wave. Typically, a borehole may be pressurized by injecting any fluid such as, for example, water, isopropyl alcohol, two-butoxyethanol, ethylene glycol, liquid nitrogen, or a slurry. It is also contemplated within the scope of the invention that a borehole may be pressurized by injecting a supercritical fluid such as, for example, supercritical carbon dioxide. Is also contemplated within the scope of the invention that a borehole may be pressurized with a gas such as, for example, atmospheric air, oxygen, nitrogen, carbon dioxide, and the like. The specific signal generator may be matched to the geological composition of the reservoir being analyzed. For example, the signal generator may be a gas in the case of a reservoir composed primarily of a gas such as natural gas or methane, while the signal generator may be a fluid or a supercritical fluid for a reservoir composed primarily of a liquid such as, for example, oil. Without being bound by theory, it is believed that matching the type of signal generator to the geological composition of the reservoir being analyzed may increase the quality and/or resolution of the data that results from the imaging of the Pf wave.

Wave propagation in the rock may be recorded with an imaging technique (such as, e.g., TFI™, SET, and the like) or a pressure monitoring technique (such as, e.g., direct measurement of pore fluid pressure with a pressure sensor deployed into a preexisting or newly drilled borehole). Poroelastic theory may then be used to estimate the permeability of the fractured rock formations from the measured wave speed. Advantageously, the present invention makes it possible to measure the permeability of fractured rock at reservoir-scale because the imaged pressure wave is able to propagate through a large rock volume, on the order of kilometers in size (e.g., 1-10 km). This is in stark contrast to traditional methods, which yield permeability measurements for much smaller rock volumes. For example, well logging tools only measure permeability in the vicinity of a borehole. Pressure transient testing accesses larger rock volumes; however, these volumes are much smaller than those enabled by the present invention, particularly in low-permeability rock formations. Another advantage of the present invention is that the signal generator may be, for example, fluid pumped into an existing borehole to allow assessment of permeability of an undisturbed reservoir (e.g., a natural reservoir that has not been subjected to artificial fracturing), or an artificially fractured reservoir (e.g., a natural reservoir that has been subjected to hydraulic fracturing, propellant fracking, or other rock stimulation methods). It is contemplated within the scope of the invention that the signal generator may be, for example, a slurry used during the process of hydraulic fracturing so that the permeability of an exemplary reservoir may be assessed at the same time that it is subject to fracturing.

The present invention provides an explanation of the physical mechanisms that give rise to the Pf wave, and demonstrates that this physical understanding may be used to estimate the permeability of the transmissive fracture network. According to the techniques herein, several field data sets show Pf waves observed via TFI™, and the Pf wave corresponds to Biot's Slow Wave. In this regard, the Biot theory may be used to estimate the permeability of the fracture networks in which the Pf waves were observed. FIG. 2 illustrates graphical views of Pf wave in the Barnett Shale. The plots suggest that the wave travels about 5,000 ft in 3 minutes, which corresponds roughly to a velocity of 6.3 m/s.

The method of the present invention includes several steps for deriving/estimating subsurface formation permeability. One step includes creation of a 3D model for seismic velocity of the rock formation under consideration. The velocity model can be constructed from sonic logs, geological models, Vertical Seismic Profiles (VSPs), crosswell tomography, seismic tomography, and the like. Another step may include installation of a surface or buried seismic array for TFI™ at the field site. Collection of passive-seismic data with the seismic array may be used to image the preexisting fracture system.

Using the 3D velocity model, the passive seismic data may be processed via TFI™ in order to create a 3D image of the seismically active fracture system. FIG. 3A illustrates a graphical view of an ambient fracture network that has been mapped by processing passive seismic data with TFI™.

A well may be drilled into the imaged rock formation, for which a permeability estimate is sought. In an exemplary embodiment, the well may be vertically oriented. Due to the attenuation of the Pf wave, it is preferable that the well be located in the center of the formation such that sufficiently high pressure perturbations can be created everywhere within the geological domain of interest.

Fluid may be injected into the borehole for some period of time at a rate such that the bottom borehole pressure increases substantially. At the same time, the rock formation is imaged with TFI™ in order to determine the fracture system that is stimulated by the generated pressure wave. While TFI™ is disclosed herein as an exemplary imaging modality, any suitable imaging system known to or conceivable by one of skill in the art could also be used. Also, the bottom borehole fluid pressure p_b is recorded as a function of time t in order to characterize the transient behavior of the source of the Pf wave. However, it is also contemplated within the scope of the invention that other variables (e.g., acceleration) could be measured in order to characterize the source of the Pf wave. FIG. 4 illustrates graphical views of borehole pressure data shows typical bottom borehole pressure data that is obtained when a slurry is injected into the well for the purposes of hydraulic fracturing.

The seismic data is processed via TFI™ in order to achieve 5D semblance data that describes the propagation of the induced pressure disturbance. The recorded 5D semblance data is used to determine the arrival time of the Pf wave front at the central surfaces of the transmissive fracture network. To that end, the data likely first need to be denoised (e.g. using anisotropic diffusion filtering), and then the front will be highlighted (e.g., using edge detection or image segmentation). FIG. 3B illustrates a graphical view of the transmissive fracture network which is illuminated by increasing the borehole pressure.

From the travel time on the central surfaces of the transmissive fracture network, the velocity of the Pf wave on the central surfaces of the transmissive fracture network is calculated. FIG. 3C illustrates a graphical view of the velocity field of the Pf wave. 5D semblance data from the passive seismic imaging is used to identify the permeable zones around the central surfaces. The rock volume in these zones contains the dominant flow paths and forms the transmissive fracture network. FIG. 3D illustrates a graphical view of the 5D semblance data from the passive seismic imaging that is used to identify the permeable zones around the major fractures.

Poroelastic theory according to Biot is used to estimate the permeability of the transmissive fracture network, k. Biot developed a theory in order to describe the propagation of elastic waves in fluid-saturated porous media. The Biot theory has been modified using the notation introduced by Stoll and Kan (R. D. Stoll and T. K. Kan. Reflection of acoustic waves at a water-sediment interface. Journal of the Acoustical Society of America, 70: 149-156, 1981) and Stoll (R. D. Stoll. Sediment Acoustics. Springer-Verlag, New York, 1989) to model parameters directly related to typical rock properties. In this regard, Biot's theory may be presented in the form of the following two partial differential equations:

$\begin{array}{cc}{\overrightarrow{\nabla }}^2\left(\mathrm{He}-C\zeta \right)=\frac{{\partial }^2}{\partial t^2}\left(\rho e-{\rho }_f\zeta \right)& \left(1\right)\\ {\overrightarrow{\nabla }}^2\left(\mathrm{Ce}-M\zeta \right)=\frac{{\partial }^2}{\partial t^2}\left({\rho }_fe-\alpha \frac{{\rho }_f}{\beta }\zeta \right)-\frac{\eta }{k}\frac{\partial \zeta }{\partial t}& \left(2\right)\end{array}$

where the volumetric strain e [−] and the increment in fluid content ξ [−] are the unknowns. The fluid pressure pf can be obtained from the primary dependent variables via

p_f=Mξ−Ce   (3)

These equations including the lumped parameters H, C and M depend on the following 11 parameters:

1. the porosity of the rock, β [−];

2. the dynamic viscosity of the fluid, η [kg·m^{31 1} s⁻¹];

3. the total mass density of the fluid-saturated porous medium, ρ [kg m⁻³];

4. the density of the fluid, ρ_f [kg m⁻³];

5. the permeability of the fractured rock, k [m²];

6. the shear modulus of the rock, μ[m²];

7. the bulk modulus of the fluid, K_f [Pa];

8. the bulk modulus of the solid material that makes up the rock that differs from the bulk modulus of the porous/fractured rock, K_r [Pa];

9. the bulk modulus of the free-draining porous/fractured rock, K_b [Pa];

10. a pore-size parameter a [m]; and

11. a factor that accounts for the tortuosity of the rock, α≧1 [−].

The lumped parameters are defined as follows:

$\begin{array}{cc}H=\frac{{\left(K_r-K_b\right)}^2}{D-K_b}+K_b+\frac{4}{3}\mu & \left(4\right)\\ C=\frac{K_r\left(K_r-K_b\right)}{D-K_b}& \left(5\right)\\ M=\frac{K_r^2}{D-K_b}& \left(6\right)\\ D=K_r\left[1+\beta \left(\frac{K_r}{K_f}-1\right)\right]& \left(7\right)\end{array}$

The following two equations of motions may be more intuitive than Equations (1) and (2):

$\begin{array}{cc}p\ddot{u}-\rho f\ddot{w}=\mu {\overrightarrow{\nabla }}^2u+\left(H-2\mu \right)\overrightarrow{\nabla }e-C\overrightarrow{\nabla }\zeta & \left(8\right)\\ C\overrightarrow{\nabla }e-M\overrightarrow{\nabla }\zeta ={\rho }_f\ddot{u}-\frac{{\rho }_f}{\beta }\ddot{w}-\frac{\eta }{k}\stackrel{.}{w}& \left(9\right)\end{array}$

where ω=β(u−U), u is the displacement of the rock skeleton, and U is the displacement of the fluid.

Equation (8) can be rewritten as follows:

(1−β)ρ_rü+βρ_fÜ=μ{right arrow over (∇)}²u+(H−2μ){right arrow over (∇)}e−C{right arrow over (∇)}ξ  (10)

where ρ_r is density of the solid material that makes up the rock. The latter equation can be interpreted as Newton's law for the rock-fluid mixture. Equation (9) can be recast as follows:

$\begin{array}{cc}{\rho }_f\ddot{U}=-\overrightarrow{\nabla }p+\frac{\eta }{k}\stackrel{.}{w}& \left(11\right)\end{array}$

Then it can be interpreted as Newton's law for the fluid where the forces per unit volume appear on the rhs of the equation, namely the negative pressure gradient and a friction term which reflects Darcy's law in case of a moving porous medium.

The problem may be assessed in Fourier space where the solutions represent harmonic waves that oscillate in time with an angular frequency ω and in space with a complex wave number 1:

e(χ, t)=B₁e^i(ωt−lχ)   (12)

ξ(χ, t)=B₂ee^i(ωt−lχ)   (13)

It is convenient to introduce a vector for the unknowns e andξ:

$\begin{array}{cc}\stackrel{->}{y}\left(x,t\right)=\left(\begin{array}{c}e\left(x,t\right)\\ \zeta \left(x,t\right)\end{array}\right)=\left(\begin{array}{c}{B}_1\\ B_2\end{array}\right)e^{\left(\mathrm{wt}-\mathrm{lx}\right)}=\stackrel{->}{B}e^{\left(\mathrm{wt}-\mathrm{lx}\right)}& \left(14\right)\end{array}$

By substituting these assumed solutions into Equations (1) and (2), one can show that the following condition must be fulfilled:

$\begin{array}{cc}\begin{array}{cc}\stackrel{\_}{H}l^2-{\mathrm{\rho \omega }}^2& {\rho }_f{\omega }^2-\stackrel{\_}{C}l^2\\ \stackrel{\_}{C}l^2-{\rho }_f{\omega }^2& m{\omega }^2-\stackrel{\_}{M}l^2-\omega F\eta /k\end{array}=0\mathrm{where}& \left(15\right)\\ F\left(\kappa \right)=\frac{1}{4}\frac{\kappa T\left(\kappa \right)}{1-\frac{2T\left(\kappa \right)}{\mathrm{\kappa }}}& \left(16\right)\end{array}$

is a function which accounts for deviations of Darcy's law in case of ω>0,

$\begin{array}{cc}T\left(\kappa \right)=\frac{{\mathrm{ber}}^{\prime }\left(\kappa \right)+{\mathrm{bei}}^{\prime }\left(\kappa \right)}{\mathrm{ber}\left(\kappa \right)+\mathrm{bei}\left(\kappa \right)}& \left(17\right)\end{array}$

is a nondimensional function, and

$\begin{array}{cc}\kappa =\alpha \sqrt{\frac{{\mathrm{\omega \rho }}_f}{\eta }}& \left(18\right)\end{array}$

is a nondimensional frequency. Equation (15) can be used to determine the dispersion relation between the complex wave number 1 and the real angular frequency ω. This dispersion relation allows one to calculate the phase velocity and the attenuation of a harmonic wave. Equation (15) can be written as a fourth-order polynomial in 1:

$\begin{array}{cc}0={\lambda }^2\underset{\underset{:=a2}{}}{\left(C^2-\mathrm{HM}\right)}+\lambda \underset{\underset{:=a1}{}}{\left(\mathrm{Hm}{\omega }^2-H\frac{\omega F\eta }{k}-2{\rho }_f{\omega }^2C+\rho M{\omega }^2\right)}+\underset{\underset{:=a0}{}}{{\rho }_f^2{\omega }^4-\rho m{\omega }^4+\mathrm{\rho }\frac{{\omega }^2F\eta }{k}}& \left(19\right)\end{array}$

where λ=12. The two solutions to this quadratic equation in λ are given by

$\begin{array}{cc}\lambda =\frac{–a_1\pm \sqrt{a_1^2-4a_0a_2}}{2a_2}& \left(20\right)\end{array}$

The two λ values correspond to Biot's first and second compressional wave. They have different phase velocities and attenuation coefficients. The complex wave number is obtained by taking the square root of each λ value. The two roots that one obtains for each λ value, l₁=l_r+il_i and l₂=−l_r−ili, represent one and the same dispersion relation. They correspond to two waves which travel with the same absolute speed into either the positive or the negative x-direction and which experience the same attenuation/growth along the direction of propagation. From Equation (15), it is also possible to calculate

$\begin{array}{cc}\stackrel{->}{B}\left(\omega \right)=\left(\begin{array}{c}{\rho }_f{\omega }^2-{\mathrm{Cl}}^2\\ {\mathrm{\rho \omega }}^2-{\mathrm{Hl}}^2\end{array}\right)& \left(21\right)\end{array}$

where l(ω) is known from the dispersion relation. Due to the linearity of the governing equations, this vector can be multiplied by any complex number. For each location x, the {right arrow over (B)} vector gives according to Equation 14 information about the ratio of the amplitudes of the e and oscillations and the phase shift between them.

From the dispersion relation, one can determine the velocity of a localized wave package, which is characterized by a dominant frequency ωc. In case of attenuation, this dominant frequency changes with time t, because the different harmonics which make up a wave package die out at different rates. In case of attenuation, the group velocity dω/dl as defined for the case without attenuation becomes complex, because 1 is complex. The real and imaginary parts of the group velocity can be related to the position of the spatial maximum of the wave package,

$\begin{array}{cc}{x}_M\left(t\right)=\mathrm{Re}{\{\frac{\omega }{l}}_{l_M}\}t& \left(22\right)\end{array}$

and the central wave number

$\begin{array}{cc}{l}_M\left(t\right)=l_c+{\Delta }^2\mathrm{Im}{\{\frac{\omega }{l}}_{l_M}\}t& \left(23\right)\end{array}$

where Δ is the width of the assumed initial Gaussian wave package. Thus the group velocity is given by

$\begin{array}{cc}{\upsilon }_g=\mathrm{Re}{\{\frac{\omega }{l}}_{l_M}\}& \left(24\right)\end{array}$

The group velocity changes with time because the central wave number 1_M changes with time according to Eq. (23).

All parameters of the Biot theory, except for the permeability k which is to be determined inversely, can be estimated for an exemplary gas shale reservoir. Shale generally has a porosity between 0 and 10%, so it is possible to estimate β=0.05. In a gas shale reservoir the fluid may typically be assumed to be methane, even though actual shale gas is a gas mixture. At 0° C. (freezing point of water) and 101.3 kPa (atmospheric pressure), methane has a viscosity η=10.2×10⁻⁶ kg m⁻¹ s⁻¹. Typical values of reservoir pressures range from 20 to 50 MPa; therefore, a gas pressure of 20 MPa may be assumed. In some instances, the reservoir pressure may be estimated by assuming a hydrostatic pressure distribution in the subsurface, in which case a pressure of 20 MPa corresponds roughly to a reservoir depth of 2 km. It may further be assumed that the reservoir temperature may be calculated from the geothermal gradient (e.g., 22.1° C. per km of depth) and that the ground surface has a temperature of 20° C., then the reservoir temperature is about 64° C. The viscosity of methane may be modeled as η=1.3×10⁻⁵ kg m⁻¹ s⁻¹. Densities of most mineral soil materials vary between 2600 and 2750 kg/m³. Even though shale is not a soil, an average value for the total density of shale, ρ=2675 kg/m³ may be used (e.g., the small porosity of the shale may be neglected in the modeling). Methane has a specific gas constant R=519 J kg⁻K⁻. Using the ideal gas law, ρf=p/(RT), ρf=114 kg/m³ is obtained for the pressure and temperature used in the viscosity calculation. A value of μ=6.1 GPa may be set according to Simm (see e.g., R. Simm Practical Gassmann fluid substitution in sand/shale sequences. First Break, 25:61-68, 2007). The bulk modulus is also the inverse of the compressibility. For an ideal gas, the isothermal compressibility is the inverse pressure. Therefore the bulk modulus is the pressure itself, which means K_f=20 MPa. K_r [Pa] is the bulk modulus of the solid material that makes up the rock that differs from the bulk modulus of the porous/fractured rock. Note that the bulk modulus is also the inverse of the compressibility. For quartz and most minerals, the compressibility is about 2×10⁻¹¹ Pa⁻¹, accordingly, a value of K_r=5×10¹⁰ Pa=50 GPa may be assumed. According to Jizba (see e.g., D. L. Jizba. Mechanical and acoustical properties of sandstones and shales. PhD thesis, Stanford University, 1991.), the bulk modulus of room-dry shale cores varies between 15 and 35 GPa. Accordingly, the arithmetic mean K_b=25 GPa may be used. In an exemplary embodiment, a=1 mm may be used for the pore-size parameter. According to Latiefa and Fauzia (F. D. E. Latiefa and U. Fauzia. Permeability and tortuosity anisotropy of a fractured porous rock. pages 1-8, 2015. 3rd International Workshop on Rock Physics, Apr. 13-17, 2015 in Perth, Western Australia.), who analyzed a 3D image of a cm-scale igneous rock sample, the tortuosity of fractured rock is on the order of 1.5. Accordingly, a value of α=1.5 may be used in an exemplary embodiment.

It is generally possible to obtain independent estimates for all Biot theory parameters of a fluid-saturated fractured rock, except for the permeability k [m²]. These parameters are the following: the porosity of the rock [−], the dynamic viscosity of the fluid [kg m⁻¹ s⁻¹], the total mass density of the fluid-saturated porous medium [kg m⁻³], the density of the fluid [kg m⁻³], the shear modulus of the rock [Pa], the bulk modulus of the fluid [Pa], the bulk modulus of the solid material that makes up the rock [Pa], the bulk modulus of the free-draining porous/fractured rock [Pa], a pore-size parameter [m], and a factor that accounts for the tortuosity of the rock [−]. Biot theory is used to calculate the dispersion relation for Biot's slow wave, i.e., the dependence of the complex wave number k on ω.

The dominant frequency ω_c of the signal generator is determined by Fourier analysis of the bottom hole pressure p_b(t), or from other variables that characterize the generator/source strength. Using the dispersion relation one can then calculate the corresponding wave number 1_c. The velocity of the Pf wave is estimated by the group velocity for a package of harmonic slow waves that is characterized by a

$\begin{array}{cc}{c}_g=r\mathrm{eal}{\{\frac{\omega }{k}}_{{\omega }_c}\}& \left(25\right)\end{array}$

Numerical values for this group velocity cannot be determined because it depends on the unknown permeability k. However, the function that describes the dependence of c_g on k can be numerically determined. FIG. 5 illustrates a graphical view of the dependence of the group velocity on permeability.

The speed of the Pf wave, c_g, in the 3D rock volume of the transmissive fracture network is used to calculate the permeability k in the same rock volume. This is done by determining at each point k(c_g) where the latter function is obtained by numerically inverting c_g(k) from Equation (25). FIG. 3E illustrates a graphical view of the calculated permeability field.

It is understood that one or more of the below methods, or aspects thereof, may be executed by at least one control unit. The term “control unit” may refer to a hardware device that includes a memory and a processor. The memory is configured to store program instructions, and the processor is specifically programmed to execute the program instructions to perform one or more processes which are described further below. Moreover, it is understood that the below methods may be executed by an apparatus comprising the control unit in conjunction with one or more other components, as would be appreciated by a person of ordinary skill in the art.

Furthermore, the control unit of the present disclosure may be embodied as non-transitory computer readable media containing executable program instructions executed by a processor, controller or the like. Examples of the computer readable mediums include, but are not limited to, ROM, RAM, compact disc (CD)-ROMs, magnetic tapes, floppy disks, flash drives, smart cards and optical data storage devices. The computer readable recording medium can also be distributed throughout a computer network so that the program instructions are stored and executed in a distributed fashion, e.g., by a telematics server or a Controller Area Network (CAN).

The steps and analysis of the present invention can be carried out using a computer, non-transitory computer readable medium, or alternately a computing device or non-transitory computer readable medium incorporated into an imaging device. Indeed, any suitable method of calculation known to or conceivable by one of skill in the art could be used. It should also be noted that while specific equations are detailed herein, variations on these equations can also be derived, and this application includes any such equation known to or conceivable by one of skill in the art. A non-transitory computer readable medium is understood to mean any article of manufacture that can be read by a computer. Such non-transitory computer readable media includes, but is not limited to, magnetic media, such as a floppy disk, flexible disk, hard disk, reel-to-reel tape, cartridge tape, cassette tape or cards, optical media such as CD-ROM, writable compact disc, magneto-optical media in disc, tape or card form, and paper media, such as punched cards and paper tape. The computing device can be a special computer designed specifically for this purpose. The computing device can be unique to the present invention and designed specifically to carry out the method of the present invention. The computing device can also take the form of an operating console computer for the imaging device. The operating console is a non-generic computer specifically designed by the scanner manufacturer for bilateral (input output) communication with the scanner. It is not a standard business or personal computer that can be purchased at a local store. Additionally this console computer carries out communications with the scanner through the execution of proprietary custom built software that is designed and written by the scanner manufacturer for the computer hardware to specifically operate the scanner hardware.

The many features and advantages of the invention are apparent from the detailed specification, and thus, it is intended by the appended claims to cover all such features and advantages of the invention, which fall within the true spirit and scope of the invention. Further, since numerous modifications and variations will readily occur to those skilled in the art, it is not desired to limit the invention to the exact construction and operation illustrated and described, and accordingly, all suitable modifications and equivalents may be resorted to, falling within the scope of the invention.

Claims

1. A method of estimating permeability of a rock formation comprising: generating, by a signal generator, a pressure (Pf) wave in a perforated portion of a borehole within the rock formation;imaging, by an imaging device, the Pf wave to produce imaging data;analyzing, by a controller, the imaging data to identify a front position of the Pf wave as a function of time and to formulate a transmissive fracture network in which the Pf wave propagates;determining the speed of the Pf wave as a function of three dimensional (3D) position within the rock formation based on the front position data; andestablishing an estimated 3D permeability field from the 3D field for the Pf wave speed using poroelastic theory.

2. The method of claim 1, wherein the signal generator is borehole pressurization.

3. The method of claim 2, wherein the borehole pressurization results from injection of a liquid, a gas, or a supercritical fluid into the borehole.

4. The method of claim 3, wherein the liquid is selected from the group consisting of water, isopropyl alcohol, two-butoxyethanol, ethylene glycol, and a slurry.

5. The method of claim 1, wherein the imaging device is a Tomographic Fracture Imaging (“TFI™”) device or a Seismic Emission Tomography (SET) device.

6. The method of claim 5, wherein the TFI™ device comprises a surface array or a sub-surface array.

7. The method of claim 6, wherein the surface array is a beam forming surface array.

8. The method of claim 6, wherein the sub-surface array is positioned 1-10 meters below the surface.

9. The method of claim 1, wherein the poroelastic theory is Biot's theory.