Fast methods to model the response of fluid-filled fractures and estimate the fracture properties

Estimation of fracture orientation and properties has become an important part of seismic reservoir characterization especially in unconventional reservoirs because of the crucial role of fractures in enhancing the permeability in tight reservoirs. The presence of fluid inside the fractures affects their seismic response. Using equivalent medium theories, seismic wave signatures such as Amplitude Variation with Offset and azimuth (AVOz), Normal Moveout (NMO) correction and shear waves splitting have been used to detect the presence of gas-filled and fluid-filled fractures. These methods, however, are unable to specify the type of fluid inside the fractures and cannot be used for thin beds and complex geology where the subsurface properties change laterally. Hence, modeling the seismic waveform using numerical methods is inevitable. The main limitation of those methods is their high computation costs. In this dissertation, I focus on developing two fast numerical methods to model the response of fluid-filled fractures as well as one fast global optimization method to estimate the fracture properties. Although local optimization methods are computationally cheap, the probability of being trapped in a local minimum becomes high when the initial model is not close to the global minimum especially when applied to highly nonlinear problems. Quantum Annealing (QA) is a recent global optimization method that was shown to be faster than Simulated Annealing (SA) in many situations. QA has been recently applied to geophysical problems. In this research, I modify QA by proposing using a new kinetic term that helps QA converge faster to the global minimum. With a synthetic dataset, I illustrate that QA is faster than Very Fast Simulated Annealing (VFSA) using a highly non-linear forward model that computes the response of seismic Amplitude Variation with Angle (AVA) for spherical waves. Most AVA inversion algorithms are based on plane wave solutions whereas seismic surveys use point sources to generate spherical waves. Although the plane wave solution is an excellent approximation for spherical waves, this approximation breaks down in the vicinity of the critical angle. Here, I implement an AVA inversion method for three parameters (P-wave velocity, S-wave velocity and density) based on analytical approximation for spherical waves. In addition, I apply this algorithm to a 2D seismic dataset from Cana field, Oklahoma with the primary objective of resolving the Woodford formation. I compare the results with those obtained by a local optimization method. The results clearly demonstrate superior performance of the proposed inversion method over that of local optimization. Specifically, the inverted images show clear delineation of the Woodford formation. For a reservoir containing vertical and rotationally invariant fractures, the linear slip model characterizes the reservoir using four properties: two elastic properties describing the isotropic host rock and two fracture properties – normal ΔN and tangential ΔT fracture weaknesses. This model, however, ignores the pore porosity effect on the anisotropy and hence the fracture properties might be inaccurate. In this work, I estimate the fracture properties as well as pore porosity using a new expression for the stiffness tensor for a porous fractured medium. I use the ray-Born approximation to calculate the seismic response of a laterally varying porous reservoir and QA to estimate the fracture properties. Using numerical experiments, I compare the inversion results from both unconstrained and constrained simultaneous (PP and PSV components) seismic inversion as well as constrained inversion using only the PP component. I explain the importance of including a constraint to mitigate the effect of the equivalence problem between ΔN and porosity. Unlike the unconstrained inversion, the estimated properties from the constrained inversion are acceptable. Also, I illustrate that the simultaneous constrained inversion is more robust than using the PP component alone. I apply this algorithm to a 3D multicomponent seismic dataset acquired in Saudi Arabia. The estimated fracture orientation agrees with those obtained in previous studies using borehole image logs, oriented cores, drilling observation and seismic in the same area. Also, the computed porosity using available well logs matches the inverted porosity very well. Computationally cheap analytical methods and equivalent medium theories available to model seismic wavefields diffracted by multiple fluid-filled fractures are not capable of handling complex fracture models or wave multi-scattering. Hence, using expensive numerical methods is inevitable. The advantages of boundary element method (BEM) over domain methods, such as finite difference and finite element methods, include the ease of handling irregular fracture geometry and reduction of the problem dimensions making the computation fast. Moreover, BEM models the complete wavefield including multiples, reverberations and refracted waves inside the fractures. The downside of BEM is that the computation cost increases rapidly whenever we increase the number of boundary elements making these methods computationally inefficient to model a large number of 2D cracks or 3D fractures. By combining the Indirect Boundary Element Method (IBEM) and a Generalized Born Series (GBS), I propose a new algorithm that can compute the response of 3D fluid-filled fracture sets effectively. In addition, when I consider equally spaced fractures that have the same geometry within a fracture set, computation can be performed even more rapidly. I compare the wavefield obtained using this approximation in five numerical experiments with those obtained from IBEM and show that the results are accurate in many situations.