CHARACTERIZING FRACTURE ORIENTATIONS IN ORTHORHOMBIC ADJACENT LAYERS USING THE PHASE OF AZIMUTHAL FOURIER COEFFICIENTS

Abstract

Methods and apparatuses characterize fracture orientations in orthorhombic adjacent layer. Seismic data with azimuthal coverage enables calculating Fourier coefficients of reflectivity at an interface between the orthorhombic adjacent layers. The phases of 2nd and 4th FCs may be used to infer the fracture orientations in the orthorhombic adjacent layers. Analysis of 2nd and 4th Fourier coefficients' phases for different incidence angles may indicate that the fracture orientations in the orthorhombic adjacent layers are aligned, orthogonal, at 45°, that one of the layers is isotropic, etc.

Description

BACKGROUND

Technical Field

Embodiments of the subject matter disclosed herein generally relate to methods and systems for seismic exploration, and, in particular, to methods and systems for investigating fracture orientations in orthorhombic media/layers.

Discussion of the Background

Seismic exploration is frequently used in oil and gas exploration and recovery for locating hydrocarbon reservoirs and designing exploitation (drilling) plans. Analysis of seismic data obtained during seismic surveys provides valuable information about structure and physical parameters of subterranean rock formations. In a seismic survey, a source generates a seismic signal that propagates into the explored formation and is, at least partially, reflected at interfaces between a formation's layer having different acoustic impedances. The source typically generates compressional seismic signals, that is, pressure variations along the propagation direction (i.e., P-waves), but both P-waves and S-waves (i.e., shear/secondary waves) may be generated at layer interfaces. Seismic detectors record seismic data which may represent both P- and S-waves.

The recorded seismic data includes a convolution of the source wavelet (i.e., the signal's signature) and the formation's response function. Seismic data processing aims to extract the formation's response function, and various parameters related to the explored formation. One way of expressing the response function is by using reflectivity, R, which for detected P-waves is R_pp. Reflectivity depends on the incidence angle, θ, (which is defined in a vertical plane) and the azimuth angle, ϕ (which is defined in a horizontal plane). Reflectivity characterizes subterranean properties and can serve as input to other inversion algorithms. For example, reflectivity can provide useful information regarding the anisotropy parameters of the subterranean formation, such as subsurface fractures. These anisotropy parameters may reveal shale plays, tight gas sands and coal bed methane, as well as carbonates in naturally fractured reservoirs. These characteristics and parameters are used to design and manage underground transportation systems, foundations of major structures, cavities for storage of liquids, gases or solids, etc.

In oil and gas exploration, this kind of information is used for determining optimal locations and orientations of wells (which may be vertical, inclined or deviated, and horizontal), to minimize wellbore instability and formation break-out. This information may also be used for stimulating the production of hydrocarbons by applying hydraulic pressure on the formation from the wellbore. Also, in geothermal exploration, such information is used to better determine location/orientation while tapping into geothermal resources. Another area where such information is usedful is in carbon capture and storage. The data inversions yield estimates of the elastic stiffnesses (velocities and anisotropic parameters) that can be used to predict lithology, porosity and the fluid content of the subsurface, as well as intensity and orientation of fractures in subterranean formations. This kind of information is important when designing drilling paths, particularly when fracking is used.

Anisotropy in orthorhombic layers (or media) is frequently due to vertical fractures parallel to one direction and, thus, perpendicular to a symmetry axis and horizontal fine layering as illustrated in FIG. 1. Adjacent layers (e.g., 110 and 120) generally have different symmetry axes, e.g., φ⁽¹⁾ and φ⁽²⁾ drawn on respective layer's top surface in FIG. 1.

It is desirable to determine the fractures' orientation from seismic data in order to better design drilling paths and/or determine locations of fracking stages.

SUMMARY

Methods and devices according to various embodiments determine fracture orientation from seismic data using the phases of azimuthal Fourier coefficients (FCs).

According to an embodiment, there is a seismic exploration method characterizing fracture orientations in orthorhombic adjacent layers. The method includes obtaining seismic data with azimuthal coverage enabling calculating azimuthal FCs of reflectivity at an interface between the orthorhombic adjacent layers, calculating phases of 2^nd and 4^th azimuthal FCs for different incidence angles, and characterizing the fracture orientations in the orthorhombic adjacent layers based on the phases of 2^nd and 4^th azimuthal FCs, for a drilling and/or fracking plan.

According to another embodiment, there is a seismic data processing apparatus including an interface and a data processing unit connected to the interface. The interface is configured to obtain seismic data with azimuthal coverage enabling calculating azimuthal FCs of reflectivity at an interface between the orthorhombic adjacent layers. The data processing unit is configured to calculate phases of 2^nd and 4^th azimuthal FCs for different incidence angles, and to characterize the fracture orientations in the orthorhombic adjacent layers based on the phases of 2^nd and 4^th azimuthal FCs, for a drilling and/or fracking plan.

According to yet another embodiment, there is a computer-readable recording medium non-transitorily storing executable codes which, when executed by a computer, make the computer perform a seismic exploration method. The method includes obtaining seismic data with azimuthal coverage enabling calculating azimuthal FCs of reflectivity at an interface between the orthorhombic adjacent layers, calculating phases of 2^nd and 4^th azimuthal FCs for different incidence angles, and characterizing the fracture orientations in the orthorhombic adjacent layers based on the phases of 2^nd and 4^th azimuthal FCs, for a drilling and/or fracking plan.

BRIEF DESCRIPTION OF THE DRAWINGS

The accompanying drawings, which are incorporated in and constitute a part of the specification, illustrate one or more embodiments and, together with the description, explain these embodiments. In the drawings:

FIG. 1 shows two adjacent orthorhombic layers with different symmetry axes of the vertical fractures;

FIG. 2 is a graph of reflectivity at an interface between two orthorhombic layers with fractures having different symmetry axes as a function of azimuth in a horizontal plane;

FIG. 3 illustrates Fourier components of reflectivity for a fixed incidence angle;

FIG. 4 illustrates the steps in relating phases of FCs to fracture orientations;

FIGS. 5 and 6 are graphical solutions for fracture orientations;

FIG. 7 is a flowchart of a method according to an embodiment;

FIG. 8 illustrates amplitude-versus-azimuth data usable to calculate FCs;

FIG. 9 is a graph illustrating phases of the 2^nd FC and of the 4^th FC as a function of the incidence angle for two orthorhombic layers not having their respective fractures aligned;

FIG. 10 is a graph illustrating phases of the 2^nd FC and of the 4^th FC as a function of the incidence angle for two orthorhombic layers with their respective fractures aligned;

FIG. 11 illustrates harmonic addition for the phase of the 2^nd FC;

FIG. 12 illustrates harmonic addition for the phase of the 4^th FC;

FIG. 13 is a graphic representation of phases; and

FIG. 14 is a schematic diagram of a computing device according to an embodiment.

DETAILED DESCRIPTION

The following description of the exemplary embodiments refers to the accompanying drawings. The same reference numbers in different drawings identify the same or similar elements. The following detailed description does not limit the invention. Reference throughout the specification to “one embodiment” or “an embodiment” means that a particular feature, structure or characteristic described in connection with an embodiment is included in at least one embodiment of the subject matter disclosed. Thus, the appearance of the phrases “in one embodiment” or “in an embodiment” in various places is not necessarily referring to the same embodiment. Further, the particular features, structures or characteristics may be combined in any suitable manner in one or more embodiments.

Embodiments described hereinafter obtain fracture orientation characterization using seismic data. More specifically, this approach uses the phase of azimuthal FCs to characterize a change of fracture orientation between two adjacent orthorhombic layers (each with substantial parallel vertical fractures embedded in a transverse isotropic media, but their parallel fractures having different symmetry axes).

Considering two adjacent orthorhombic layers as illustrated in FIG. 1, the reflectivity R_pp varies with azimuth ϕ at this interface as illustrated in FIG. 2, exhibiting peaks and troughs. Reflectivity may be decomposed in a Fourier series

R(Φ,θ)=Σ_n=0^N(u_n(θ)cos(nΦ)+v_n(θ)sin(nΦ)).  (1)

that can be rewritten and simplified by retaining only dominant terms as

R(Φ,θ)=r₀(θ)+r₂(θ)cos 2(Φ−Φ₂(θ))r₄(θ)cos 4)(Φ−Φ₄(θ))  (2)

where θ is the incidence angle in the vertical plane at the reflection/refraction point of the interface and Φ is the azimuth angle in a horizontal plane at the reflection/refraction point of the interface. FIG. 3 illustrates Fourier components of reflectivity R₀, R₂ and R₄ for a fixed value of the incident angle, that is: R₀=r₀; R₂ (Φ)=r₂ cos 2(Φ−Φ₂); and R₄(Φ)=r₄ cos 4(Φ−Φ₄).

Evaluating anisotropy has lately been included in most workflows for seismic processing and reservoir characterization. The gain over workflows considering media isotropic is undeniable but the increased complexity associated with taking into consideration anisotropy remains a challenge. Lowering the symmetry of the anisotropy from transverse isotropic (TI) to orthorhombic to monoclinic implies increasing the number of parameters. Therefore, simplifying the problem by reducing the dimensionality of the parameterization is important.

A powerful tool to simplify the complexity of anisotropic reflection coefficients are azimuthal Fourier coefficients, FC, (described, e.g., in “Amplitude variations with azimuths (AVAZ) inversion based on Linearized inversion of common azimuthal sections” by L. Ikelle, included in “Seismic Anisotropy” edited by E. Fjr et al., SEG 1996, pp. 601-644; the 2001 article “Azimuth-dependent AVO in reservoirs containing non-orthogonal fracture sets” by Sayers, C. and Dean, S. published in Geophysical Prospecting, No. 49, pp. 100-106; and the article “Azimuthal Fourier Coefficients” by Downton, J., Roure, B. and Hunt, L. published in CSEG Reporter of December 2011, pp. 22-36, all of which are incorporated herewith by reference in their entirety). In the 2015 article “Interpreting azimuthal FCs for anisotropic and fracture parameters” by Downton, J. and Roure, B. published in Interpretation of August 2015, pp. ST9-ST27 (incorporated herewith by reference in its entirety), it is shown that azimuthal FCs allow decoupling anisotropic reflection coefficients into sub-problems associated with different types of anisotropy, hence reducing cross-talk between parameters.

Sayers and Dean′ article (already specified) interprets the phase of FC for monoclinic anisotropy in terms of linear slip theory thereby enabling to relate changes of phase to differences of compliance ratio between the various fracture sets. The issue is relating phases of FCs to fracture orientations and vice versa.

FIG. 4 illustrates the steps in establishing the theoretical basis of relating phases of FCs to fracture orientations. At 410, the orthorhombic stiffness matrices are defined for both top and bottom layers. As shown in the 1997 article “Orthorhombic media: Modeling elastic wave behavior in a vertically fractured earth” by Schoenberg, M. and Helbig, K., published in Geophysics, Vol. 62, No. 6, pp. 1954-1974 (which is incorporated herewith by reference in its entirety), a density-normalized stiffness matrix for a single set of parallel vertical fractures perpendicular to the x-axis in a background transverse isotropic (TI) media is given by

$\begin{array}{cc}a=\left(\begin{array}{cccccc}{a}_{11}& a_{12}& a_{13}& & & \\ & a_{22}& a_{23}& & & \\ & & a_{33}& & & \\ & & & a_{44}& & \\ & & & & a_{55}& \\ & & & & & a_{66}\end{array}\right)& \left(4\right)\end{array}$

which matrix relates the six independent components of the strain tensor to the independent components of the stress tensor.

Initially, the horizontal symmetry axis is chosen along the x-axis. In order to consider the actual orientation of the symmetry axis for each layer (i.e., ϕ⁽¹⁾ for the top and ϕ⁽²⁾ for the bottom layer as illustrated in FIG. 1), each density-normalized stiffness matrix is rotated relative to its symmetry axis using a Bond transformation at step 420. Bond transformation, which is described in the 1990 article “Velocity anisotropy terminology for geophysicists” by Winterstein, D. F., published in Geophysics, Vol. 55, No. 8, pp. 1070-1088 (incorporated herewith by reference in its entirety), is a counter-clockwise rotation about the z-axis yielding 6×6 matrices a′ for each layer with the following non-zero elements (in the following formula ϕ is ϕ⁽¹⁾ or ϕ⁽²⁾ as appropriate):

a′ ₃₃ =a ₃₃

a′ ₁₁=(a₆₆+½a₁₂)sin²2Φ+a₁₁ cos⁴Φ+a₂₂ sin⁴Φ

a′ ₁₂=¼(−4a₆₆+a₁₁+a₂₂)sin²2Φ+a₁₂(cos⁴Φ+sin⁴Φ)

a′ ₂₂=(a₆₆+½a₁₂)sin²2Φ+a₂₂ cos⁴Φ+a₁₁ sin⁴Φ

a′ ₁₆=−½a₆₆ sin 4Φ+½a₁₁ sin 2Φ cos²Φ−½a₂₂ sin 2Φ sin²Φ−½a₁₂ sin²2Φ

a′ ₂₆=½a₆₆ sin 4Φ+½a₁₁ sin 2Φ sin²Φ−½a₂₂ sin 2Φ cos²Φ+½a₁₂ sin²2Φ

a′ ₆₆ =a ₆₆ cos²2Φ+¼(a₁₁+a₂₂−2a₁₂)sin²2Φ

a′ ₁₃ =a ₁₃ cos²Φ+a₂₃ sin²Φ

a′ ₂₃ =a ₂₃ cos²Φ+a₁₃ sin²Φ

a′ ₄₄ =a ₄₄ cos²Φ+a₅₅ sin²Φ

a′ ₅₅ =a ₅₅ cos²Φ+a₄₄ sin²Φ

a′ ₃₆=½(a₁₃−a₂₃)sin 2Φ

a′ ₄₅=½(a₅₅−a₄₄)sin 2Φ  (5)

Then, at 430, the weak anisotropy parameters and their contrasts for general anisotropy are calculated in the manner introduced in the 1998 article “Polarization, phase velocity, and NMO velocity of qP-waves in arbitrary weakly anisotropic media” by Pšenčik, I. and Gajewski, D., published in Geophysics, Vol. 63, No. 5, pp. 1754-1766, and the 2001 article “Properties of weak contrast PP reflection/transmission coefficients for weakly anisotropic elastic media” by Pšenčik, I. and Martins, J. L., published in Studia Geophysica et Geodaetica, No. 45, pp. 176-199 (which are incorporated herewith by reference in their entirety). The weak anisotropy parameters for both layers are:

$\hspace{1em}\begin{array}{cc}\{\begin{array}{c}{\varepsilon }_x=\frac{a_{11}^{\prime }-{\alpha }^2}{2{\alpha }^2},{\varepsilon }_y=\frac{a_{22}^{\prime }-{\alpha }^2}{2{\alpha }^2},{\varepsilon }_z=\frac{a_{33}^{\prime }-{\alpha }^2}{2{\alpha }^2}\\ {\delta }_x=\frac{a_{13}^{\prime }+2a_{55}^{\prime }-{\alpha }^2}{{\alpha }^2},{\delta }_y=\frac{a_{23}^{\prime }+2a_{44}^{\prime }-{\alpha }^2}{{\alpha }^2},{\delta }_z=\frac{a_{12}^{\prime }+2a_{66}^{\prime }-{\alpha }^2}{{\alpha }^2}\\ {\chi }_x=\frac{a_{14}^{\prime }+2a_{56}^{\prime }}{{\alpha }^2},{\chi }_y=\frac{a_{25}^{\prime }+2a_{46}^{\prime }}{{\alpha }^2},{\chi }_z=\frac{a_{36}^{\prime }+2a_{45}^{\prime }}{{\alpha }^2}\\ {\varepsilon }_{15}=\frac{a_{15}^{\prime }}{{\alpha }^2},{\varepsilon }_{16}=\frac{a_{16}^{\prime }}{{\alpha }^2},{\varepsilon }_{24}=\frac{a_{24}^{\prime }}{{\alpha }^2},{\varepsilon }_{26}=\frac{a_{26}^{\prime }}{{\alpha }^2},{\varepsilon }_{34}=\frac{a_{34}^{\prime }}{{\alpha }^2},{\varepsilon }_{35}=\frac{a_{35}^{\prime }}{{\alpha }^2}\\ {\gamma }_x=\frac{a_{55}^{\prime }-{\beta }^2}{2{\beta }^2},{\gamma }_y=\frac{a_{44}^{\prime }-{\beta }^2}{2{\beta }^2},{\varepsilon }_{45}=\frac{a_{45}^{\prime }}{{\alpha }^2}\end{array}& \left(6\right)\end{array}$

where α and β are the reference P-wave and S-wave velocities, respectively. The parameter contrasts are then:

$\hspace{1em}\begin{array}{cc}\{\begin{array}{c}{\mathrm{\Delta \varepsilon }}_x={\varepsilon }_x^{\left(2\right)}-{\varepsilon }_x^{\left(1\right)},\Delta {\varepsilon }_y={\varepsilon }_y^{\left(2\right)}-{\varepsilon }_y^{\left(1\right)},\Delta {\varepsilon }_z={\varepsilon }_z^{\left(2\right)}-{\varepsilon }_z^{\left(1\right)}\\ {\mathrm{\Delta \delta }}_x={\delta }_x^{\left(2\right)}-{\delta }_x^{\left(1\right)},\Delta {\delta }_y={\delta }_y^{\left(2\right)}-{\delta }_y^{\left(1\right)},\Delta {\delta }_z={\delta }_z^{\left(2\right)}-{\delta }_z^{\left(1\right)}\\ {\mathrm{\Delta \gamma }}_x={\gamma }_x^{\left(2\right)}-{\gamma }_x^{\left(1\right)},\Delta {\gamma }_y={\gamma }_y^{\left(2\right)}-{\gamma }_y^{\left(1\right)},\Delta {\chi }_z={\chi }_z^{\left(2\right)}-{\chi }_z^{\left(1\right)}\\ \Delta {\varepsilon }_{16}={\varepsilon }_{16}^{\left(2\right)}-{\varepsilon }_{16}^{\left(1\right)},\Delta {\varepsilon }_{26}={\varepsilon }_{26}^{\left(2\right)}-{\varepsilon }_{26}^{\left(1\right)},\Delta {\varepsilon }_{45}={\varepsilon }_{45}^{\left(2\right)}-{\varepsilon }_{45}^{\left(1\right)}\end{array}& \left(7\right)\end{array}$

where the upper scripts (1) and (2) refer to the top and bottom layers, respectively.

Further, at 440, FCs are calculated using the weak anisotropy contrasts as follows (see formulas in the 2015 article by Downton):

$\begin{array}{cc}\{\begin{array}{c}{B}_{C2}=\frac{1}{4}\left(\Delta {\delta }_x-\Delta {\delta }_y\right)-2\stackrel{\_}{g}\left(\Delta {\gamma }_x-\Delta {\gamma }_y\right)\\ B_{S2}=\frac{1}{2}\left({\mathrm{\Delta \chi }}_z-4\stackrel{\_}{g}\Delta {\varepsilon }_{45}\right)\\ C_{C2}=\frac{1}{4}\left(\Delta {\varepsilon }_x-\Delta {\varepsilon }_y\right)\\ C_{S2}=\frac{1}{4}\left(\Delta {\varepsilon }_{16}+\Delta {\varepsilon }_{26}\right)\\ C_{C4}=\frac{1}{16}\left(\Delta {\varepsilon }_x+\Delta {\varepsilon }_y-\Delta {\delta }_z\right)\\ C_{S4}=\frac{1}{8}\left(\Delta {\varepsilon }_{16}-\Delta {\varepsilon }_{26}\right)\end{array}& \left(8\right)\end{array}$

where g=β²/α² and the overbar indicates average of the properties across the interface. The FCs for a given incidence angle θ are then:

$\begin{array}{cc}\{\begin{array}{c}{u}_2=\left(B_{C2}+C_{C2}{\tan }^2\theta \right){\sin }^2\theta \\ v_2=\left(B_{S2}+C_{S2}{\tan }^2\theta \right){\sin }^2\theta \\ u_4=C_{C4}{\sin }^2\theta {\tan }^2\theta \\ v_4=C_{S4}{\sin }^2\theta {\tan }^2\theta \end{array}& \left(9\right)\end{array}$

the FCs phases are then defined as:

$\begin{array}{cc}\{\begin{array}{c}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{v_2}{u_2}\right)\\ {\varphi }_4=\frac{1}{4}\arctan \left(\frac{v_4}{u_4}\right)\end{array}.& \left(10\right)\end{array}$

At 450, the phases of the 2^nd and the 4^th FCs are expressed as function or the symmetry axes angles ϕ⁽¹⁾ and ϕ⁽²⁾. By replacing formula (9) into (10), the following equations are obtained for the phases of the 2^nd and the 4^th FCs:

$\begin{array}{cc}\{\begin{array}{c}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{\left(B_{S2}+{\tan }^2\theta C_{S2}\right){\sin }^2\theta }{\left(B_{C2}+{\tan }^2\theta C_{C2}\right){\sin }^2\theta }\right)=\frac{1}{2}\arctan \left(\frac{C_{S2}}{C_{C2}}\left(\frac{B_{S2}\text{/}{C}_{S2}+{\tan }^2\theta }{B_{C2}\text{/}{C}_{C2}+{\tan }^2\theta }\right)\right)\\ {\varphi }_4=\frac{1}{4}\arctan \left(\frac{C_{S4}{\sin }^2\theta {\tan }^2\theta }{C_{C4}{\sin }^2\theta {\tan }^2\theta }\right)=\frac{1}{4}\arctan \left(\frac{C_{S4}}{C_{C4}}\right)\end{array}.& \left(11\right)\end{array}$

The phase of the 2^nd FC varies with the incidence angle, whereas the phase of the 4^th FC does not. In the special case where B_S2/C_S2=B_C2/C_C2, the 2^nd FC loses its incidence angle dependency too. If formulas (5) to (11) are combined, the phases of the 2^nd and the 4^th FCs can be expressed directly as functions of the orthorhombic media parameters and their horizontal symmetry axes angles ϕ⁽ⁱ⁾:

$\begin{array}{cc}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{\left(X^{\left(2\right)}+Y^{\left(2\right)}{\tan }^2\theta \right)\sin 2{\varphi }^{\left(2\right)}-\left(X^{\left(1\right)}+Y^{\left(1\right)}{\tan }^2\theta \right)\sin 2{\varphi }^{\left(1\right)}}{\left(X^{\left(2\right)}+Y^{\left(2\right)}{\tan }^2\theta \right)\cos 2{\varphi }^{\left(2\right)}-\left(X^{\left(1\right)}+Y^{\left(1\right)}{\tan }^2\theta \right)\cos 2{\varphi }^{\left(1\right)}}\right)\mathrm{and}& \left(12\right)\\ {\varphi }_4=\frac{1}{4}\arctan \left(\frac{Z^{\left(2\right)}\sin 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\sin 4{\varphi }^{\left(1\right)}}{Z^{\left(2\right)}\cos 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\cos 4{\varphi }^{\left(1\right)}}\right).\mathrm{where}& \left(13\right)\\ \{\begin{array}{c}X\left(i\right)=\frac{1}{4}\frac{a_{13}^{\left(i\right)}-a_{23}^{\left(i\right)}}{{\alpha }^{\left(i\right)2}}+\frac{1}{2}\left(a_{55}^{\left(i\right)}-a_{44}^{\left(i\right)}\right)\left(\frac{1}{{\alpha }^{\left(i\right)2}}-\frac{2g}{{\beta }^{\left(i\right)2}}\right)\\ Y^{\left(i\right)}=\frac{1}{8}\frac{a_{11}^{\left(i\right)}-a_{22}^{\left(i\right)}}{{\alpha }^{\left(i\right)2}}\\ Z\left(i\right)=\frac{1}{16}\frac{1}{{\alpha }^{\left(i\right)2}}\left(\frac{1}{2}\left(a_{11}^{\left(i\right)}+a_{22}^{\left(i\right)}\right)-a_{12}^{\left(i\right)}-2a_{66}^{\left(i\right)}\right).\end{array}& \left(14\right)\end{array}$

For orthorhombic media, X⁽ⁱ⁾, Y⁽ⁱ⁾ and Z⁽ⁱ⁾ are typically non-zero. The phase of the 2^nd FC does not depend on the incidence angle if equations (16) using notations (15) are satisfied:

$\begin{array}{cc}\{\begin{array}{c}{B}_{C2}=X^{\left(2\right)}\cos 2{\varphi }^{\left(2\right)}-X^{\left(1\right)}\cos 2{\varphi }^{\left(1\right)}\\ B_{S2}=X^{\left(2\right)}\sin 2{\varphi }^{\left(2\right)}-X^{\left(1\right)}\sin 2{\varphi }^{\left(1\right)}\\ C_{C2}=Y^{\left(2\right)}\cos 2{\varphi }^{\left(2\right)}-Y^{\left(1\right)}\cos 2{\varphi }^{\left(1\right)}\\ C_{S2}=Y^{\left(2\right)}\sin 2{\varphi }^{\left(2\right)}-Y^{\left(1\right)}\sin 2{\varphi }^{\left(1\right)}\\ C_{C4}=Z^{\left(2\right)}\cos 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\cos 4{\varphi }^{\left(1\right)}\\ C_{S4}=Z^{\left(2\right)}\sin 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\sin 4{\varphi }^{\left(1\right)}\end{array}& \left(15\right)\\ \frac{B_{S2}}{C_{S2}}=\frac{B_{C2}}{C_{C2}}\iff \frac{X^{\left(2\right)}\sin 2{\varphi }^{\left(2\right)}-X^{\left(1\right)}\sin 2{\varphi }^{\left(1\right)}}{Y^{\left(2\right)}\sin 2{\varphi }^{\left(2\right)}-Y^{\left(1\right)}\sin 2{\varphi }^{\left(1\right)}}=\frac{X^{\left(2\right)}\cos 2{\varphi }^{\left(2\right)}-X^{\left(1\right)}\cos 2{\varphi }^{\left(1\right)}}{Y^{\left(2\right)}\cos 2{\varphi }^{\left(2\right)}-Y^{\left(1\right)}\cos 2{\varphi }^{\left(1\right)}}.& \left(16\right)\end{array}$

The symmetry axes ϕ⁽¹⁾ and ϕ⁽²⁾ that solve equation (16) are sought for known values of X⁽ⁱ⁾, Y⁽ⁱ⁾ and Z⁽ⁱ⁾ (i=1, 2). FIGS. 5 and 6 represent equation (16) in nuances of gray with ϕ⁽¹⁾ and ϕ⁽²⁾ along the axes. FIG. 5 is a linear scale and FIG. 6 a logarithmic scale making particular solutions more visible, such as when fractures on top and bottom layers are aligned and when they are orthogonal.

Thus, if the fractures on the top and the bottom layers are aligned, that is, ϕ⁽¹⁾=ϕ⁽²⁾ then

$\begin{array}{cc}\{\begin{array}{c}\frac{B_{S2}}{C_{S2}}=\frac{X^{\left(2\right)}-X^{\left(1\right)}}{Y^{\left(2\right)}-Y^{\left(1\right)}}\\ \frac{B_{C2}}{C_{C2}}=\frac{X^{\left(2\right)}-X^{\left(1\right)}}{Y^{\left(2\right)}-Y^{\left(1\right)}}.\end{array}& \left(17\right)\end{array}$

If the fractures on the top layer are orthogonal on the ones in the bottom layer, i.e., ϕ⁽¹⁾=ϕ⁽²⁾+90°, then

$\begin{array}{cc}\{\begin{array}{c}\frac{B_{S2}}{C_{S2}}=\frac{X^{\left(2\right)}+X^{\left(1\right)}}{Y^{\left(2\right)}+Y^{\left(1\right)}}\\ \frac{B_{C2}}{C_{C2}}=\frac{X^{\left(2\right)}+X^{\left(1\right)}}{Y^{\left(2\right)}+Y^{\left(1\right)}}.\end{array}& \left(18\right)\end{array}$

Both in the case of aligned and in the case of orthogonal fractures, equation (11) becomes

$\begin{array}{cc}\{\begin{array}{c}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{C_{S2}}{C_{C2}}\right)\\ {\varphi }_4=\frac{1}{4}\arctan \left(\frac{C_{S4}}{C_{C4}}\right).\end{array}& \left(19\right)\end{array}$

If C_S2=C_C2=0, then the phase of the 2^nd FC no longer depends on the incidence angle, that is:

$\begin{array}{cc}\{\begin{array}{c}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{B_{S2}}{B_{C2}}\right)\\ {\varphi }_4=\frac{1}{4}\arctan \left(\frac{C_{S4}}{C_{C4}}\right).\end{array}& \left(20\right)\end{array}$

If, besides C_S2=C_C2=0, the fractures are aligned or orthogonal, then ϕ₂=ϕ₄=ϕ⁽²⁾.

In view of these formulas, FIG. 7 is a flowchart of a method 700 for characterizing fracture orientations in orthorhombic adjacent layers according to an embodiment. Method 700 includes obtaining seismic data with azimuthal coverage enabling calculating FCs at 710. For example, seismic data may be azimuthal angle stacks or any other processed land or marine data allowing the computation of FCs. FIG. 8 illustrates amplitude-versus-azimuth data usable to calculate FCs for incidence angles 0°, 10°, 20°, 30° and 40°; that is, curves 800, 810, 820, 830 and 840, respectively. In FIG. 8, the circles mark amplitude magnitudes of 0.1, 0.2, etc. and data azimuth 0.1) (e.g., between source and receiver) varies counterclockwise as marked on the circumference.

Thus, the FCs and their phases are calculated at 720. FIGS. 9 and 10 are graphs illustrating phases of the 2^nd FC (the continuous line) and of the 4^th FC as functions of the incidence angle. The fractures in the two orthorhombic media for which phases of the FCs are shown in FIG. 9 are not aligned, while the fractures in the two orthorhombic media for which phases of the FCs are shown in FIG. 10 are aligned. Thus, when the fractures are aligned, the phase of the 2^nd FC no longer depends on the incident angle and it is equal to the phase of the 4^th FC and the fracture orientation in both layers.

Returning now to FIG. 7, at 730, the fracture orientations in the orthorhombic adjacent layers are characterized based on the phases of the 2^nd and 4^th FCs. This characterization can be approached from different points of view.

One way to approach this characterization is by using a signal interpretation of the phases. Based on the harmonic addition theorem, the phase of an FC is the phase of the sum of two signals whose respective phases are the symmetry axes of the two orthorhombic layers. The phase of the 2^nd FC is the phase of a signal S₂ that is a sum of signals S₂⁽¹⁾ and S₂⁽²⁾ as follows:

$\begin{array}{cc}\{\begin{array}{c}{S}_2^{\left(1\right)}\left(\varphi \right)=-\left(X^{\left(1\right)}+Y^{\left(1\right)}{\tan }^2\theta \right)\sin 2\left(\varphi -{\varphi }^{\left(1\right)}\right)\\ S_2^{\left(2\right)}\left(\varphi \right)=\left(X^{\left(2\right)}+Y^{\left(2\right)}{\tan }^2\theta \right)\sin 2\left(\varphi -{\varphi }^{\left(2\right)}\right)\\ S_2\left(\varphi \right)=S_2^{\left(1\right)}\left(\varphi \right)+S_2^{\left(2\right)}\left(\varphi \right)=Y\sin 2\left(\varphi -{\varphi }_2\right)\end{array}& \left(21\right)\end{array}$

where

Y ²=(X⁽¹⁾+Y⁽¹⁾ tan²θ)²+(X⁽²⁾+Y⁽²⁾ tan²θ)²−2(X⁽¹⁾+Y⁽¹⁾ tan²θ) . . . (X⁽²⁾+Y⁽²⁾ tan²θ)cos 2(ϕ⁽²⁾−ϕ⁽¹⁾).

FIG. 11 illustrates these signals as functions of azimuth ϕ. The intersections of S₂⁽¹⁾, S₂ and S₂⁽²⁾ with the x axis yield ϕ⁽¹⁾, ϕ₂ and ϕ⁽²⁾, respectively.

Similarly, the phase of the 4^th FC is the phase of a signal S₄ that is a sum of signals S₄⁽¹⁾ and S₄⁽²⁾ as follows:

$\begin{array}{cc}\{\begin{array}{c}{S}_4^{\left(1\right)}\left(\varphi \right)=-Z^{\left(1\right)}\sin 4\left(\varphi -{\varphi }^{\left(1\right)}\right)\\ S_4^{\left(2\right)}\left(\varphi \right)=Z^{\left(2\right)}\sin 4\left(\varphi -{\varphi }^{\left(2\right)}\right)\\ S_4\left(\varphi \right)=S_4^{\left(1\right)}\left(\varphi \right)+S_4^{\left(2\right)}\left(\varphi \right)=Z\sin 4\left(\varphi -{\varphi }_4\right)\end{array}& \left(22\right)\end{array}$

where

Z ² =Z ⁽¹⁾² +Z ⁽²⁾²−2Z⁽¹⁾Z⁽²⁾ cos 4(ϕ⁽²⁾−ϕ⁽¹⁾).

FIG. 12 illustrates S₄⁽¹⁾, S₄ and S₄⁽²⁾ as functions of azimuth ϕ. The intersections of S₄⁽¹⁾, S₄ and S₄⁽²⁾ with the x axis yield ϕ⁽¹⁾, ϕ₄ and ϕ⁽²⁾, respectively.

Another way to relate FC phases to symmetry axes is graphically illustrated in FIG. 13. The phase of the Fourier coefficients is geometrically related to the symmetry axes, scaled by the periodicity (or order) of the Fourier coefficients and the anisotropy of each orthorhombic medium. Formulas (12) and (13) can be generically re-written as

$\begin{array}{cc}\tan \left(n{\varphi }_n\right)=\frac{x^{\left(2\right)}\sin \left(n{\varphi }^{\left(2\right)}\right)-x^{\left(1\right)}\sin \left(n{\varphi }^{\left(1\right)}\right)}{x^{\left(2\right)}\cos \left(n{\varphi }^{\left(2\right)}\right)-x^{\left(1\right)}\cos \left(n{\varphi }^{\left(1\right)}\right)}& \left(23\right)\end{array}$

with n=2 or 4. Thus, on a generic circle, one starts by marking angles nϕ⁽¹⁾ and nϕ⁽²⁾. The horizontal and vertical projections of nϕ⁽¹⁾ and nϕ⁽²⁾ are sin(nϕ⁽¹⁾), cos(nϕ⁽¹⁾), sin(nϕ⁽²⁾) and cos(nϕ⁽²⁾), respectively. Then, on the y axis are drawn points corresponding to x⁽¹⁾ sin(nϕ⁽¹⁾) and to x⁽²⁾ sin(nϕ⁽²⁾) and on the x axis are drawn points corresponding to x⁽¹⁾ cos(nϕ⁽¹⁾) and to x⁽²⁾ cos(nϕ⁽²⁾). Point A of coordinates (x⁽²⁾ cos(nϕ⁽²⁾),x⁽²⁾ sin(nϕ⁽²⁾)) is at angle nϕ⁽²⁾ and point B of coordinates (x⁽¹⁾ cos(nϕ⁽¹⁾), x⁽¹⁾ sin(nϕ⁽¹⁾)) is at angle nϕ⁽¹⁾. The difference between the projection of point A on the y axis and the projection of point B on the y axis is equal to segment CA and represents the numerator in formula (23). The difference between the projection of point A on the x axis and the projection of point B on the x axis is equal to segment B₁C=−BC and represents the denominator in formula (23). Thus, angle AB₁C is nϕ_n its tan function being equal to formula (23).

If the symmetry axes angles and the phases of FCs are small enough (e.g., less than 30° not valid in many cases), then a first order Taylor approximation of (23) yields:

$\begin{array}{cc}n{\varphi }_n=\frac{x^{\left(2\right)}n{\varphi }^{\left(2\right)}-x^{\left(1\right)}n{\varphi }^{\left(1\right)}}{x^{\left(2\right)}-x^{\left(1\right)}}& \end{array}$

which is equivalent with

$\begin{array}{cc}{\varphi }_n=\frac{x^{\left(2\right)}{\varphi }^{\left(2\right)}-x^{\left(1\right)}{\varphi }^{\left(1\right)}}{x^{\left(2\right)}-x^{\left(1\right)}}.& \left(24\right)\end{array}$

Thus, in this case, the phases of the FCs are weighted averages of the symmetry axes angles, weighted by the anisotropy in each layer.

In the case one layer is isotropic and the other is orthorhombic, the incidence angle dependency of the phase of the 2^nd FC is controlled by the ratios B_S2/C_S2 and B_C2/C_C2:

$\begin{array}{cc}\{\begin{array}{c}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{\left(B_{S2}+{\tan }^2\theta C_{S2}\right){\sin }^2\theta }{\left(B_{C2}+{\tan }^2\theta C_{C2}\right){\sin }^2\theta }\right)=\frac{1}{2}\arctan \left(\frac{C_{S2}}{C_{C2}}\left(\frac{B_{S2}/C_{S2}+{\tan }^2\theta }{B_{C2}/C_{C2}+{\tan }^2\theta }\right)\right)\\ {\varphi }_4=\frac{1}{4}\arctan \left(\frac{C_{S4}{\sin }^2\theta {\tan }^2\theta }{C_{C4}{\sin }^2\theta {\tan }^2\theta }\right)=\frac{1}{4}\arctan \left(\frac{C_{S4}}{C_{C4}}\right).\end{array}& \left(25\right)\end{array}$

If one of the top or bottom layer is isotropic and the other orthorhombic, then B_S2/C_S2=B_C2/C_C2=X⁽²⁾/Y⁽²⁾ and the phase of the 2^nd FC does not depend on the incident angle, the phases of the 2^nd and 4^th FCs being equal to the orthorhombic symmetry axis angle.

If there are the same properties in both layers except for the fracture orientations, then the phase of the 2^nd FC does not depend on the incident angle. Using the properties of the inverse trigonometric functions, the phases are the average of the symmetry axes shifted by a constant:

$\begin{array}{cc}\{\begin{array}{c}{\varphi }_2=\pm \frac{\pi }{4}+\frac{1}{2}\left({\varphi }^{\left(2\right)}+{\varphi }^{\left(1\right)}\right)\\ {\varphi }_4=\pm \frac{\pi }{8}+\frac{1}{2}\left({\varphi }^{\left(2\right)}+{\varphi }^{\left(1\right)}\right)\end{array}& \left(26\right)\end{array}$

While the phase of the 4^th FC is only sensitive to the anisotropy (anellipticity) in the far angles, the phase of the 2^nd FC is sensitive to anisotropy both in the mid and far angles, giving rise to its incident angle dependency. If there is no anisotropy in the far angles (i.e., Y(i)=0) or if the anisotropy is the same in the mid and far angles (i.e. X(i)=Y(i)), then the angle dependency of the phase of the 2^nd FC disappears:

$\begin{array}{cc}{\varphi }_2=\frac{1}{2}\arctan \left(\frac{X^{\left(2\right)}\sin 2{\varphi }^{\left(2\right)}-X^{\left(1\right)}\sin 2{\varphi }^{\left(1\right)}}{X^{\left(2\right)}\cos 2{\varphi }^{\left(2\right)}-X^{\left(1\right)}\mathrm{co}s2{\varphi }^{\left(1\right)}}\right).& \left(27\right)\end{array}$

If the fracture orientations are 45° apart, then the phase of the 2^nd FC remains dependent on incidence angle, but the phase of the 4^th FC is equal to the symmetry axes (which are the same modulo 45°):

$\begin{array}{cc}{\varphi }_4=\frac{1}{4}\arctan \left(\frac{Z^{\left(2\right)}\sin 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\sin 4{\varphi }^{\left(1\right)}}{Z^{\left(2\right)}\cos 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\cos 4{\varphi }^{\left(1\right)}}\right){\varphi }_4=\frac{1}{4}\arctan \left(\frac{Z^{\left(2\right)}\sin 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\sin 4\left({\varphi }^{\left(2\right)}+\frac{\pi }{4}\right)}{Z^{\left(2\right)}\cos 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\cos 4\left({\varphi }^{\left(2\right)}+\frac{\pi }{4}\right)}\right){\varphi }_4=\frac{1}{4}\arctan \left(\frac{Z^{\left(2\right)}\sin 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\sin \left(4{\varphi }^{\left(2\right)}+\pi \right)}{Z^{\left(2\right)}\cos 4{\varphi }^{\left(2\right)}-Z^{\left(1\right)}\cos \left(4{\varphi }^{\left(2\right)}+\pi \right)}\right){\varphi }_4=\frac{1}{4}\arctan \left(\frac{Z^{\left(2\right)}\sin 4{\varphi }^{\left(2\right)}+Z^{\left(1\right)}\sin 4{\varphi }^{\left(2\right)}}{Z^{\left(2\right)}\cos 4{\varphi }^{\left(2\right)}+Z^{\left(1\right)}\cos 4{\varphi }^{\left(2\right)}}\right){\varphi }_4=\frac{1}{4}\arctan \left(\frac{\sin 4{\varphi }^{\left(2\right)}{Z}^{\left(2\right)}+Z^{\left(1\right)}}{\cos 4{\varphi }^{\left(2\right)}{Z}^{\left(2\right)}+Z^{\left(1\right)}}\right){\varphi }_4=\frac{1}{4}\arctan \left(\tan 4{\varphi }^{\left(2\right)}\right)={\varphi }^{\left(2\right)}.& \left(28\right)\end{array}$

The above-discussed methods may be implemented in a computing device 1400 as illustrated in FIG. 14. Hardware, firmware, software or a combination thereof may be used to perform the various steps and operations described herein.

Exemplary computing device 1400 suitable for performing the activities described in the exemplary embodiments may include a server 1401. Server 1401 may include a central processor (CPU) 1402 coupled to a random-access memory (RAM) 1404 and to a read-only memory (ROM) 1406. ROM 1406 may also be other types of storage media to store programs, such as programmable ROM (PROM), erasable PROM (EPROM), etc. Processor 1402 may communicate with other internal and external components through input/output (I/O) circuitry 1408 and bussing 1410 to provide control signals and the like. Processor 1402 carries out a variety of functions as are known in the art, as dictated by software and/or firmware instructions.

Server 1401 may also include one or more data storage devices, including hard drives 1412, CD-ROM drives 1414 and other hardware capable of reading and/or storing information, such as DVD, etc. In one embodiment, software for carrying out the above-discussed steps may be stored and distributed on a CD-ROM or DVD 1416, a USB storage device 1418 or other form of media capable of portably storing information. These storage media may be inserted into, and read by, devices such as CD-ROM drive 1414, disk drive 1412, etc. Server 1401 may be coupled to a display 1420, which may be any type of known display or presentation screen, such as LCD, plasma display, cathode ray tube (CRT), etc. A user input interface 1422 is provided, including one or more user interface mechanisms such as a mouse, keyboard, microphone, touchpad, touch screen, voice-recognition system, etc.

Server 1401 may be coupled to other devices, such as sources, detectors, etc. The server may be part of a larger network configuration as in a global area network (GAN) such as the internet 1428, which allows ultimate connection to various computing devices.

According to one embodiment, I/O circuitry 1408 is configured to obtain seismic data with azimuthal coverage enabling calculating FCs of reflectivity at an interface between the orthorhombic adjacent layers, and processor 1402 is configured to calculate phases of 2^nd and 4^th azimuthal FCs for different incidence angles and to characterize the fracture orientations in the orthorhombic adjacent layers based on the phases of 2^nd and 4^th azimuthal FCs, for a drilling and/or fracking plan.

In yet another embodiment, RAM 1404 stores executable codes that, when executed, make the I/O circuitry 1408 obtain seismic data with azimuthal coverage enabling calculation of FCs of reflectivity at an interface between the orthorhombic adjacent layers, and processor 1402 to calculate phases of the 2^nd and 4^th azimuthal FCs for different incidence angles and to characterize the fracture orientations in the orthorhombic adjacent layers based on the phases of the 2^nd and 4^th azimuthal FCs, for a drilling and/or fracking plan.

The disclosed embodiments provide methods and apparatuses for characterizing fracture orientations in orthorhombic adjacent layers based on azimuthal Fourier coefficients of reflectivity at an interface between the orthorhombic adjacent layers. It should be understood that this description is not intended to limit the invention. On the contrary, the exemplary embodiments are intended to cover alternatives, modifications and equivalents, which are included in the spirit and scope of the invention as defined by the appended claims. Further, in the detailed description of the exemplary embodiments, numerous specific details are set forth in order to provide a comprehensive understanding of the claimed invention. However, one skilled in the art would understand that various embodiments may be practiced without such specific details.

Although the features and elements of the present embodiments are described in the embodiments in particular combinations, each feature or element can be used alone without the other features and elements of the embodiments or in various combinations with or without other features and elements disclosed herein. The methods or flowcharts provided in the present application may be implemented in a computer program, software or firmware tangibly embodied in a computer-readable storage medium for execution by a general-purpose computer or a processor.

This written description uses examples of the subject matter disclosed to enable any person skilled in the art to practice the same, including making and using any devices or systems and performing any incorporated methods. The patentable scope of the subject matter is defined by the claims, and may include other examples that occur to those skilled in the art. Such other examples are intended to be within the scope of the claims.

Claims

1. A seismic exploration method characterizing fracture orientations in orthorhombic adjacent layers, the method comprising: obtaining seismic data with azimuthal coverage enabling calculating azimuthal Fourier Coefficients, FCs, of reflectivity at an interface between the orthorhombic adjacent layers;calculating phases of 2nd and 4th azimuthal FCs for different incidence angles; andcharacterizing the fracture orientations in the orthorhombic adjacent layers based on the phases of 2nd and 4th azimuthal FCs, for a drilling and/or fracking plan.

2. The method of claim 1, wherein the phases of 2nd and 4th azimuthal FCs indicate whether the fracture orientations in the orthorhombic adjacent layers are aligned.

3. The method of claim 1, wherein the phases of 2nd and 4th azimuthal FCs indicate whether the fracture orientations in the orthorhombic adjacent layers are orthogonal on each other.

4. The method of claim 1, wherein azimuth angle dependences of the phases of 2nd and 4th azimuthal FCs are used to infer the fracture orientations in the orthorhombic adjacent layers.

5. The method of claim 1, wherein if one of the orthorhombic adjacent layers is an isotropic layer, the phases of 2nd and 4th azimuthal FCs are substantially equal, independent of the incidence angle and substantially equal to an angle of a symmetry axis of other among the orthorhombic adjacent layers.

6. The method of claim 1, wherein if the orthorhombic adjacent layers have same properties except the fracture orientations, then the phases of the azimuthal FCs are averages of the fracture orientations shifted by respective constants.

7. The method of claim 1, wherein, if the fracture orientations are 45° apart, the phase of the 2nd FC remains dependent on incidence angle but the phase of the 4th FC is equal to the symmetry axes which are same if expressed modulo 45°.

8. The method of claim 1, wherein the seismic data is made of azimuthal angle stacks for at least two incident angles.

9. A seismic data processing apparatus, comprising: an interface configured to obtain seismic data with azimuthal coverage enabling calculating azimuthal Fourier Coefficients, FCs, of reflectivity at an interface between the orthorhombic adjacent layers; anda data processing unit connected to the interface and configured to calculate phases of 2nd and 4th azimuthal FCs for different incidence angles; andto characterize the fracture orientations in the orthorhombic adjacent layers based on the phases of 2nd and 4th azimuthal FCs, for a drilling and/or fracking plan.

10. The seismic data processing apparatus of claim 9, wherein the phases of 2nd and 4th azimuthal FCs indicate whether the fracture orientations in the orthorhombic adjacent layers are aligned.

11. The seismic data processing apparatus of claim 9, wherein the phases of 2nd and 4th azimuthal FCs indicate whether the fracture orientations in the orthorhombic adjacent layers are orthogonal on each other.

12. The seismic data processing apparatus of claim 9, wherein azimuth angle dependence of the phases of 2nd and/or 4th azimuthal FCs are used to infer the fracture orientations in the orthorhombic adjacent layers.

13. The seismic data processing apparatus of claim 9, wherein if one of the orthorhombic adjacent layers is an isotropic layer, the phases of 2nd and 4th azimuthal FCs are substantially equal, independent of the incidence angle and substantially equal to an angle of a symmetry axis of other among the orthorhombic adjacent layers.

14. The seismic data processing apparatus of claim 9, wherein if the orthorhombic adjacent layers have same properties except the fracture orientations, then the phases of the 2nd and 4th azimuthal FCs are averages of the fracture orientations shifted by respective constants.

15. The seismic data processing apparatus of claim 9, wherein the seismic data is made of azimuthal angle stacks for at least two incident angles.

16. A computer readable recording medium non-transitorily storing executable codes that when executed by a computer make computer perform a method for characterizing fracture orientations in orthorhombic adjacent layers, the method comprising: obtaining seismic data with azimuthal coverage enabling calculating azimuthal Fourier Coefficients, FCs, of reflectivity at an interface between the orthorhombic adjacent layers;calculating phases of 2nd and 4th azimuthal FCs for different incidence angles; andcharacterizing the fracture orientations in the orthorhombic adjacent layers based on the phases of 2nd and 4th azimuthal FCs, for a drilling and/or fracking plan.

17. The computer readable recording medium of claim 16, wherein the phases of 2nd and 4th azimuthal FCs indicate whether the fracture orientations in the orthorhombic adjacent layers are aligned.

18. The computer readable recording medium of claim 16, wherein the phases of 2nd and 4th azimuthal FCs indicate whether the fracture orientations in the orthorhombic adjacent layers are orthogonal on each other.

19. The computer readable recording medium of claim 16, wherein azimuth angle dependence of the phases of 2nd and/or 4th azimuthal FCs are used to infer the fracture orientations in the orthorhombic adjacent layers.

20. The computer readable recording medium of claim 16, wherein if one of the orthorhombic adjacent layers is an isotropic layer, the phases of 2nd and 4th azimuthal FCs are substantially equal, independent of the incidence angle and substantially equal to an angle of a symmetry axis of other among the orthorhombic adjacent layers.