A practical approach to develop a proper anisotropic rock physics model for media comprising of multiple fracture sets in the absence of sophisticated laboratory/wireline data

For economical production from a fractured reservoir, a characteristic analysis of the fracture parameters like its density and orientation within the reservoir is essential to improve the fluid flow during extraction. This study deals with the development of a proper anisotropic rock physics model for a media with multiple fracture sets to study the spatial distribution of important fracture parameters i.e., fracture density and orientation in the absence of sophisticated laboratory/wireline and pre-stack seismic data. The crest of hydrocarbon producing fault-bounded Balkassar Anticline in Northern Potwar, Upper Indus Basin, Pakistan is selected as a case study representing a potential zone for development of fractures at reservoir level (Sakesar Limestone). The methodology consists of the interpretation of 3D post-stack seismic and conventional wireline log data to demarcate the reservoir containing fractures. The Ant-tracking discrete fracture network (DFN) attribute is applied on 3D post-stack seismic data to obtain an initial estimate about the presence of fracture corridors and their orientations. Based on this initial estimate, a proper rock physics model has been developed utilizing inverse Gassmann relations, T-matrix approximation, and Brown and Korringa relations. The output from the developed rock physics model has been displayed in the form of 13 effective independent elastic stiffness constants (monoclinic symmetry–representing media comprising of multiple fracture sets) as a function of fracture densities and azimuthal fracture orientations. A clear decreasing trend in effective elastic stiffness constants with increasing fracture densities can be observed. Similarly, a periodic trend of effective elastic stiffness constants with fracture orientations can be observed. These trends are more or less expected, but they would have been difficult to quantify without a proper rock physics model. The use of independent effective elastic constants for the generation of synthetic seismic amplitude versus angle and azimuth (AVAZ) data and its correlation with observed seismic AVAZ data in a geostatistical sense has been discussed.

1 Introduction

Fractures are commonly found in carbonate reservoir rocks (Zhiqian et al., 2016; Bagni et al., 2020) and hold the key for the extraction of hydrocarbons from the reservoir. Prolific management of the fractured reservoirs involves effectively defining fluid flow pathways during production for which effective characteristic analysis of fracture systems is essential (Ali and Jakobsen, 2011a; Ali and Jakobsen, 2011b). Furthermore, fracture properties such as fracture density and azimuthal fracture orientation, when studied spatially, can assist in the optimization of fractured reservoir production (Sayers, 2009; Ali et al., 2015). However, despite the advancement in methodologies pertaining to subsurface data acquisition and processing, identification and adequate characteristic analysis of fractures crucial to production remain a task yet to be accomplished (Gholipour et al., 2016; Bagni et al., 2020).

The investigation of fractures within a carbonate reservoir can be carried out efficiently through sophisticated subsurface borehole data (Bagni et al., 2020). However, the fracture corridors within the carbonate reservoirs are typically extended over tens to hundreds of meters in terms of height and width while being regionally spread over a kilometer (Huang et al., 2017). Furthermore, the well-based petrophysical analysis used in the characterization of fracture properties (like fracture density) cannot be available throughout the reservoir (Jadoon et al., 2006; Bagni et al., 2020). Therefore, characterizing a fractured reservoir using borehole data may require drilling more wells, which is practically not feasible. In order to tackle this problem, a simple yet practical approach is to develop a proper rock physics model which can be utilized to generate synthetic (calculated) seismic amplitude versus angle and azimuth (AVAZ) data. This synthetic seismic AVAZ data can be correlated with observed seismic AVAZ data, if available, in a geostatistical sense to obtain the spatial distribution of important fracture parameters required to enhance fluid flow.

Rock physics modeling is a vital tool for studying the seismic-based characterization of a reservoir (Misaghi et al., 2010). The focus of rock physics concepts is to link the composition of a rock with its macroscopic properties (Wang et al., 2018). The use of an appropriate rock physics model can, therefore, delineate characteristics of a rock effectively (Wang et al., 2018). In terms of seismic fracture characterization, rock physics modelling offers a reduced number of unknown parameters, thus offering a relative advantage during seismic inversion (Ali and Jakobsen, 2014). Furthermore, a proper rock physics model delivers a better physical description of the subsurface provided that the model is able to relate the fracture parameters with the effective elastic stiffness constants (Ali and Jakobsen, 2014).

The reason for the difficulty in the prediction of subsurface fracture systems from seismic data is the lack of resolution at which these fractures can be mapped (Souque et al., 2019). The accurate prediction of the fracture systems requires sophisticated laboratory/wireline and pre-stack seismic data, including core cuttings, full bore formation micro imager (FMI) logs, and seismic amplitude versus angle and azimuth (AVAZ) data at the reservoir level. However, this data is not easily available. Therefore, it is a challenging task to develop a proper anisotropic rock physics model utilizing conventional wireline log and 3D post-stack seismic data.

This study proposes a practical approach for the development of a proper rock physics model in the case of a reservoir (Sakesar Formation in this study) containing multiple fracture sets in the absence of sophisticated laboratory/wireline data. The crest of the Balkassar anticline has been taken as the case study as it provides a geological advantage of high probability of presence of multiple fracture sets. Balkassar oil field is present in Potwar basin, a north-western part of Himalayan Mountain ranges (Gee and Gee, 1989; Kemal et al., 1991) as depicted in Figure 1. A compressional regime acting on the area led to the development of anticlinal structural traps bounded by reverse faults. Patala Shale and Pre-Cambrian Salt Range Formation are considered to be the main source rocks in the locality (Bender et al., 1995; Masood et al., 2017). The main reservoirs are the fractured carbonates of the Sakesar and Chorgali formations in the Balkassar oil field with the Chorgali Formation also acting as the seal rock (Khan et al., 1986).

FIGURE 1

FIGURE 1. Location of Balkassar Area in Kohat-Potwar Plateau on a map depicting regional structures in detail (Gee and Gee, 1989).

The input data for developing the rock physics model is obtained from available conventional wireline log data (three tracks–lithology, resistivity, and porosity) of Well OXY-01 and 3D seismic post-stack data acquired at the Balkassar oil field. The workflow followed for the generation of such a rock physics model in the case of multiple fracture sets is displayed in Figure 2. section 2 of this paper focuses of the methodology followed to achieve the objective of developing a proper rock physics model. section 3 discusses the results obtained from each applied method including the output of the proposed rock physics model. Section 4 deals with the discussion and implication of the applied method, while section 5 concludes the study.

FIGURE 2

FIGURE 2. Workflow followed in the study for development of a rock physics model to characterize a medium with multiple fracture sets.

2 Methodology

2.1 Demarcation of fractured reservoir (Sakesar Limestone)

Interpretation of the 3D seismic line focuses on the extension of the formation as well as the demarcation of regional structures. In this study, 3D seismic data of Balkassar area with a dominant frequency of 35 Hz has been interpreted for the purpose of developing the discrete fracture network (DFN) model at the crest of Balkassar anticline (Figure 3). The important horizons including Chorgali, Sakesar, Patala, and Lockhart formations have been marked. The demarcation has been done with the help of synthetic seismogram (seismic to well tie tool) generated from the well log data of OXY-01 and wavelet extracted from the seismic traces nearest to the well. Faults have been marked on the basis of seismic discontinuities and background geological knowledge (Figure 3).

FIGURE 3

FIGURE 3. Interpretation of 3D seismic Inline 235 showing Balkassar Anticline bound by reverse faults. Formations have been marked on the basis of synthetic seismogram tied with OXY-01 well.

2.2 Discrete fracture network modeling

Model development of complete 3D characteristic profiling of fracture networks present within naturally fractured reservoirs (NFR) is essential for improving the supervision of a reservoir and consequently, the rate of recovery (Narr et al., 2006; Bisdom et al., 2014). DFN modeling is a proficient tool used for the depiction of the spatial distribution of fractures as well as the details of properties like location, size, density, orientation, and conductivity of individual fractures within an NFR (Tran et al., 2006). It takes into account the geometry of fractures, their conductive ability, and interconnection to construct a fracture network (Aydin and Akin, 2019). Furthermore, DFN ensures the contribution of each fracture as a separate entity unlike the continuum approach and dual-porosity modeling. It, therefore, establishes a more realistic display of fracture systems and paves a path to study response to the flow impact from properties of individual fractures (Aydin and Akin, 2019). This study makes use of the DFN model for Sakesar Formation in the Balkassar area to ensure the type of fractures present in the subsurface reservoir before proceeding to develop a rock physics model based on this priori information.

The Ant-Tracking algorithm has been utilized for the construction of the DFN model. The ant-tracking algorithm simulates a fracture pattern analogous to ants making pathways optimal for searching food (Hu et al., 2014). During their search for food, ants leave pheromones along the pathways to inform other ants about the possibly shortest passage to any located food (Silva et al., 2005). The same concept is used in the Ant-tracking algorithm for 3D seismic volume (Pedersen et al., 2002). Ant-tracking distributes a vast number of electronic “ants” over a seismic volume with each “ant” moving through a fractured surface-emitting “pheromone” for the rest to follow (Hu et al., 2014). Hence a large number of “ants” will trace planar surfaces like faults and fractures in a 3D seismic volume marked by the “pheromones” that the “ants” left compared to the non-structure element i.e., noise, which will be marked by a smaller number of “pheromones” (Cox and Seitz, 2007; Hu et al., 2014). Furthermore, deviation up to only 15° is allowed to “ants” from their original direction in order to improve and reinforce the tracing of planar structure (faults and fractures) rather than seismic irregularity due to noise (Fang et al., 2017).

Ant-tracking process uses a series of geophysical attributes which are inserted sequentially in a background program for the detection of faults and fractures on a 3D seismic volume (Ouenes, 2000). The general workflow to run the Ant-track volume is given in Figure 4. The final output is an attribute volume that shows faults and fracture zones in detail. In this workflow, amplitude contrast with dip guide has been used as an edge detection attribute. Two attributes have been computed in series from the result of the edge detection i.e., 3D edge enhancement and structural smoothing (Figure 5). Output from edge detection is subtracted from the input seismic data. The edge evidence attribute is computed on the output of edge enhancement which further enhances the discontinuities of interest. The final fault volume is obtained by multiplying the output of edge enhancement (A), with the output of edge evidence (B), by using the following syntax in the seismic calculator: if (A>0, A, 0)B. On Fault Cube, Ant track attribute slice is generated to extract regional faults visible on seismic data. The results obtained from the Ant track attribute applied to the seismic data are given in Figures 5A,B.

FIGURE 4

FIGURE 4. Generalized workflow followed for Ant-Tracking DFN model generation used for identification of fracture network system.

FIGURE 5

FIGURE 5. (A) Fracture analysis of Sakesar Formation using Ant-tracking attribute developed for 3D seismic volume at Balkassar anticline (B) The time-slice of the Ant-tracking attribute at Sakesar level depicting the presence of multiple fracture sets.

The application of DFN Ant-track algorithm on 3D seismic post-stack data at the crest of Balkassar anticline gives us an initial estimate of fracture distribution. This initial estimate helps in defining the conceptual model for fracture distribution at the reservoir level, within the Sakesar Formation, for which the effective anisotropic elastic properties can be obtained utilizing rock physics algorithms.

2.3 Rock physics modeling

In carbonates rocks, fractures play imperative role in providing a suitable pathway for flow of the hydrocarbons. The primary/secondary migration and alignment of these fractures to preferred orientation lead us towards the anisotropic behavior of subsurface rock formation (Ali and Jakobsen, 2011a; Ali and Jakobsen, 2011b). Therefore, in order to study the characteristics of fractures in fractured media, an anisotropic rock physics model has to be developed.

Multiple fracture sets in Sakesar Formation have been identified by means of a seismic Ant-tracking algorithm. On the basis of this observation, a monoclinic anisotropic rock physics model, taking fracture density and fracture orientation as the variables under observation has been developed. The workflow for rock physics modeling is shown in Figure 6. The monoclinic anisotropic system gives 13 independent effective elastic stiffness constants which have been plotted against the fracture densities and fracture orientations of the fractures sets. These effective elastic constants are required for the computation of AVAZ data at the reservoir level.

FIGURE 6

FIGURE 6. Scheme followed in this research for developing a rock physics model in case of multiple fracture sets.

2.3.1 Mineral properties

The Sakesar Formation acting as a reservoir in the study area consists of compact limestone as dominant lithology indicated by the petrophysical analysis of Well OXY-01, with shale volume of 4%, effective porosity of 5%, and water saturation of 27.7% (Table 1; Figure 7). The standard values of bulk and shear modulus of mineral (calcite) along with density is taken from (Mavko et al., 2009) listed in Table 1.

TABLE 1

TABLE 1. Input data for rock physics algorithms obtained from petrophysical analysis of OXY-01 well at reservoir (Sakesar Formation) level (Figure 7). The bulk and shear modulus of mineral (calcite) along with density is taken from (Mavko et al., 2009). The in-situ velocities and density obtained from wireline log data helps to obtain saturated bulk modulus ($K_{sat}$) at reservoir level for computation of dry rock properties.

FIGURE 7

FIGURE 7. Petrophysical analysis of the Well OXY-01 depicting the zone of interest and the overburden zone within the compact Sakesar Formation.

2.3.2 Dry rock matrix properties (without fractures)

The bulk modulus for the frame ($K_{dry}$) of the reservoir can be measured from laboratory experimentation, empirical relation, or formulation using wireline log data (Kumar, 2006). $K_{dry}$ can be determined from wireline log data using rephrased Gassmann equation for $K_{dry}$ (Zhu and McMechan, 1990; Smith et al., 2003; Kumar, 2006; Mavko et al., 2009) given in Eq. 1 as:

$K_{dry}=\frac{K_{sat}\left(\frac{\varnothing K_{mat}}{K_{fl}}+1-\varnothing \right)-K_{mat}}{\frac{\varnothing K_{mat}}{K_{fl}}+\frac{K_{sat}}{K_{mat}}-1-\varnothing }.(1)$

Here, $K_{dry}$ is bulk modulus of dry porous rock, $K_{sat}$ is the insitu value of the bulk modulus of saturated rock, $K_{mat}$ is matrix (mineral) bulk modulus, $K_{fl}$ is effective fluid bulk modulus and $ø$ is porosity. The $K_{sat}$ is found using the density ($\rho$), P-wave ($V_p$), and S-wave ($V_S$) logs using the relation given in Eq. 2 as:

$K_{sat}=V_P^{\mathit{2}}\ast \rho -\frac{\mathit{4}}{\mathit{3}}{\mu }_{\mathrm{s}\mathrm{a}\mathrm{t}}(2)$

where, ${\mu }_{\mathrm{s}\mathrm{a}\mathrm{t}}$ is the insitu shear modulus which is equal to ${\rho V}_S^{\mathit{2}}$. $V_P$ and $V_S$ (in m/s) are estimated by taking inverse of compressional slowness and shear slowness respectively. $K_{fl}$ is calculated for a fractured media saturated with oil and water using Wood’s equation for homogenous saturation (Wood, 1955) given in Eq. 3 as:

$\frac{\mathit{1}}{K_{fl}}=\frac{S_w}{K_w}+\frac{S_o}{Ko}(3)$

where, $K_w$ and$K_o$ are bulk moduli of water and oil, respectively, whereas $S_w$ and $S_o$ are saturation of water and oil respectively. The $S_w$ has been estimated using resistivity and porosity logs of Well OXY-01. The $S_o$ can be estimated by subtracting the $S_w$ from 1, since the total sum of saturation of the two fluids is equal to 1. In particular, the bulk moduli for fluid phases (water and oil in this case) have been computed using (Batzle and Wang, 1992) relations. The detailed workflow for computing the bulk moduli of fluid phases is provided by Kumar (2006). The values used as input for the determination of $K_{fl}$and $K_{dry}$ are given in Table 1.

2.3.3 T-matrix for incorporating the effect of fractures in poroelastic monoclinic model

For the dry case of the fractured porous medium, the dry effective stiffness tensor ${\mathit{C}}_{\mathit{d}}^{*}$ has been estimated using the T-matrix approach shown in Eq. 4, given by (Jakobsen et al., 2003a; Jakobsen et al., 2003b; Ali and Jakobsen, 2011a; Ali and Jakobsen, 2011b; Ali et al., 2015).

${\mathit{C}}_{\mathit{d}}^{\ast }={\mathit{C}}^{\left(\mathit{0}\right)}+{\mathit{C}}_{\mathit{1}}:{\left({\mathit{I}}_{\mathit{4}}+{\mathit{C}}_{\mathit{1}}^{-\mathit{1}}:{\mathit{C}}_{\mathit{2}}\right)}^{-\mathit{1}},(4)$

where, ‘$:$’ implies double scalar product (Auld, 1990), ${\mathit{C}}^{\left(\mathit{0}\right)}$ adds the matrix properties input into the equation and is termed as the background stiffness tensor, ${\mathit{I}}_{\mathit{4}}$ is the identity matrix for the fourth-rank tensors, and ${\mathit{C}}_{\mathit{1}}$ is the fourth rank tensor representing first-order corrections in order to incorporate the effects of isolated fractures. ${\mathit{C}}_{\mathit{2}}$ is the second-order correction which incorporates the fracture-fracture interaction effect. ${\mathit{C}}_{\mathit{1}}$ and ${\mathit{C}}_{\mathit{2}}$ are computed using Eqs 5, 7.

${\mathit{C}}_{\mathit{1}}={{\displaystyle \sum }}_{r=\mathit{1}}^N{v}^{\left(r\right)}{\mathit{t}}_{\mathit{d}}^{\left(\mathit{r}\right)},(5)$

where, $v^{\left(r\right)}$ is the porosity of inclusion, ${\mathit{t}}_{\mathit{d}}^{\left(\mathit{r}\right)}$ incorporates the effects of individual fractures of elastic stiffness ${\mathit{C}}^{\left(\mathit{0}\right)}$ and is given in Eq. 6 as:

${\mathit{t}}_{\mathit{d}}^{\left(\mathit{r}\right)}={-\mathit{C}}^{\left(\mathit{0}\right)}:{\left[{\mathit{I}}_{\mathit{4}}+{\mathit{G}}^{\left(\mathit{r}\right)}:{\mathit{C}}^{\left(\mathit{0}\right)}\right]}^{-\mathit{1}},(6)$

where, ${\mathit{G}}^{\left(\mathit{r}\right)}$ is a fourth-rank tensor integrated over a characteristic spheroid with shape similar to fractures of type $r$ presented by strain green’s function (Jakobsen et al., 2003a; Jakobsen et al., 2003b; Ali and Jakobsen, 2011a; Ali and Jakobsen, 2011b; Ali et al., 2015). The second-order correction is given by:

${\mathit{C}}_{\mathit{2}}={\displaystyle \sum }_{r=\mathit{1}}^N{{\displaystyle \sum }}_{s=\mathit{1}}^N{v}^{\left(r\right)}{\mathit{t}}_{\mathit{d}}^{\left(r\right)}:{\mathit{G}}_{\mathit{d}}^{\left(rs\right)}:{\mathit{t}}_{\mathit{d}}^{\left(s\right)}{v}^{\left(s\right)}.(7)$

Here, ${\mathit{G}}_{\mathit{d}}^{\left(\mathit{r}\mathit{s}\right)}$ is the strain green’s function in the form of a fourth-rank tensor which represents the spatial distribution of fractures over an ellipsoid determining the symmetry of the correlation for fractures (Carcione, 1995; Ali et al., 2015). This function estimates the probability of distribution of the two fracture sets namely $r$ and $s$. The function ${\mathit{G}}_{\mathit{d}}^{\left(\mathit{r}\mathit{s}\right)}$ is integrated over an ellipsoid with an aspect ratio equal to ${\mathrm{p}}^{\left(\mathrm{s}|\mathrm{r}\right)}\left(\mathrm{x}-{\mathrm{x}}^{\mathrm{ʹ}}\right)$, that calculates the probability density for determining a type $\mathrm{s}$ inclusion at ${\mathrm{x}}^{\mathrm{ʹ}}$ provided a type r inclusion is present at the point $\mathrm{x}$ (Jakobsen et al., 2003a; Jakobsen et al., 2003b; Ali and Jakobsen, 2011a; Ali and Jakobsen, 2011b; Ali et al., 2015). The correlation function ${\mathrm{p}}^{\left(\mathrm{s}|\mathrm{r}\right)}\left(\mathrm{x}-{\mathrm{x}}^{\mathrm{ʹ}}\right)$ defines how individual fractures are distributed throughout the strata (Nguyen and Nam, 2011; Ali et al., 2015). The symmetry of the correlation function is represented by the aspect ratio. A rough aspect ratio value used in this study was measured to be 1/1,000 from exposed strata of Sakesar Formation in the vicinity of acquired seismic data. The porosity $v^{\left(r\right)}$ of type $r$ fractures relate to the fracture density ${\epsilon }^{\left(r\right)}$ through the equation $v^{\left(r\right)}=\frac{\mathit{4}}{\mathit{3}}{\epsilon }^{\left(r\right)}{\alpha }^{\left(r\right)}$ where ${\alpha }^{\left(r\right)}$ represents the type $r$ fracture’s aspect ratio. The modified form of ${\mathit{t}}_{\mathit{d}}^{\left(\mathit{r}\right)}$ for the two sets of fracture parameters can now be written as given in Eq. 8.

${\mathit{t}}_{\mathit{d}}^{\left(\mathit{r}\right)}={\mathit{t}}_{\mathit{d}}^{\left(\mathit{r}\right)}\left({\psi }_{\mathit{1}},{\psi }_{\mathit{2}},{\epsilon }_{\mathit{1}},{\epsilon }_{\mathit{2}}\right),(8)$

where, ${\psi }_{\mathit{1}}$ , ${\psi }_{\mathit{2}}$ denotes the azimuthal fracture orientations for the two fracture sets while ${\epsilon }_{\mathit{1}}$ , ${\epsilon }_{\mathit{2}}$ are the fracture densities for the respective fracture set. In general, T-matrix gives the initial output in form of VTI (Vertically transversely isotropic) model. This output matrix upon 90° rotation gives HTI (Horizontally transversely isotropic) symmetry. An arbitrary rotation through bond transformation converts the output matrix into monoclinic symmetry representing the multiple fracture sets (Figure 6). The effect of fluid in this anisotropic poroelastic model has been incorporated using an anisotropic form of Gassmann’s equation called Brown-Korringa relation (Brown and Korringa, 1975) which can be written in the symbolic matrix notation given in Eq. 9 as:

${\mathit{S}}^{*}={\mathit{S}}_{\mathit{d}}^{*}+\frac{\left({\mathit{S}}_{\mathit{d}}^{*}-{\mathit{S}}_{\mathit{m}}\right):\left({\mathit{I}}_{\mathbf{2}}{\otimes \mathit{I}}_{\mathbf{2}}\right):\left({\mathit{S}}_{\mathit{d}}^{*}-{\mathit{S}}_{\mathit{m}}\right)}{{\varphi }^0\left({\mathit{I}}_{\mathbf{2}}.{\mathit{S}}_{\mathit{m}}.{\mathit{I}}_{\mathbf{2}}-\frac{1}{K_{fl}}\right)-{\mathit{I}}_{\mathbf{2}}.\left({\mathit{S}}_{\mathit{d}}^{*}-{\mathit{S}}_{\mathit{m}}\right).{\mathit{I}}_{\mathbf{2}}}(9)$

Here $\otimes$ symbolizes the dyadic product of tensor (Auld, 1990), ${\mathit{I}}_{\mathbf{2}}$ is second rank tensor identity, ${\mathit{S}}_{\mathit{m}}$ represents the compliance tensor for solid mineral constituent (calcite), ${\mathit{S}}_{\mathit{d}}^{*}$ denotes the effective compliance tensor for dry fractured rock and ${\mathit{S}}^{*}$ inculcates the effective compliance tensor for the saturated fractured porous medium into the equation. Total porosity is represented by ${\varphi }^0$ which is the sum of storage porosity of the homogenous matrix and the fracture porosity.

The final product of rock physics modeling for two mesoscopic fracture sets with different fracture orientations and fracture densities for monoclinic symmetry is 13 independent effective elastic stiffness constants. These 13 independent constants in terms of Voigt index notation are ${\mathit{C}}_{\mathit{11}},{\mathit{C}}_{\mathit{22}},{\mathit{C}}_{\mathit{33}},{\mathit{C}}_{\mathit{12}},{\mathit{C}}_{\mathit{13}},{\mathit{C}}_{\mathit{23}},{\mathit{C}}_{\mathit{44}},{\mathit{C}}_{\mathit{55}},{\mathit{C}}_{\mathit{66}},{\mathit{C}}_{\mathit{16}},{\mathit{C}}_{\mathit{26}},{\mathit{C}}_{\mathit{36}}$ _, and ${\mathit{C}}_{\mathit{45}}$_. The 13 independent constants have been varied with fracture orientations of the two fracture sets while keeping their fracture densities constant and vice versa. The Voigt condensed notation in the matrix from for the case of monoclinic assumption having 13 independent stiffness constants is given in Eq. 10 as:

$C_{ij}=\left[\begin{array}{cccccc}{C}_{\mathit{11}}& C_{\mathit{12}}& C_{\mathit{13}}& \mathit{0}& C_{\mathit{15}}& \mathit{0}\\ C_{\mathit{12}}& C_{\mathit{22}}& C_{\mathit{23}}& \mathit{0}& C_{\mathit{25}}& \mathit{0}\\ C_{\mathit{13}}& C_{\mathit{23}}& C_{\mathit{33}}& \mathit{0}& C_{\mathit{35}}& \mathit{0}\\ \mathit{0}& \mathit{0}& \mathit{0}& C_{\mathit{44}}& \mathit{0}& C_{\mathit{46}}\\ C_{\mathit{15}}& C_{\mathit{25}}& C_{\mathit{35}}& \mathit{0}& C_{\mathit{55}}& \mathit{0}\\ \mathit{0}& \mathit{0}& \mathit{0}& C_{\mathit{46}}& \mathit{0}& C_{\mathit{66}}\end{array}\right](10)$

3 Results

3.1 Spatial extension of Sakesar Formation

A synthetic seismogram is generated using sonic and density logs of OXY-01 for the purpose of seismic-to-well-tie. With the help of this, the top and bottom of Sakesar Formation have been marked at 1.37 and 1.475 s respectively, by correlating the seismic (time domain) data to well (depth domain) data. The interpretation of seismic Inline 235 (Figure 3) shows a clear anticlinal structure which is called the Balkassar anticline. Reverse faults formed as a result of slippage caused by compressional regime bounding the anticline have also been interpreted (Figure 3).

3.2 Ant-tracking DFN model analysis

The DFN model effectively delineates the fractures present within the carbonate reservoir. Greater values displayed by Ant-tracking volume shows high fractured zone while low values correspond to the less fractured strata which can be related directly with secondary porosity evolution in carbonates (Figure 5A). Fractures are concentrated more towards the south within the anticline while in general, the highly populated fractures towards the east corresponds to the faulting phenomena towards the east of the anticline (Figure 5A). Two dominant fracture sets, almost orthogonal to one another, oriented in NE-SW and NW-SE direction can be observed in the formation (Figure 5B).

3.3 Rock physics model for monoclinic symmetry

A rock physics model for monoclinic symmetry has been developed for Sakesar Limestone reservoir in Balkassar locality. 13 elastic stiffness constants have been computed and displayed against fracture orientations (${\psi }_{\mathit{1}}$ ,${\psi }_{\mathit{2}}$) and fracture densities (${\epsilon }_{\mathit{1}},{\epsilon }_{\mathit{2}}$) of two mesoscopic fracture sets observed in the formation using Ant-tracking DFN model (Figure 5B). Elastic stiffness constants as a function of ${\epsilon }_{\mathit{1}}$ and ${\epsilon }_{\mathit{2}}$ have been plotted in Figure 8 keeping ${\psi }_{\mathit{1}}$ and ${\psi }_{\mathit{2}}$ constant at 35° and 45° respectively. From Figure 8, it is clear that there is a decrease in elastic stiffness constants with increasing fracture densities showing an inverse linear relation.

FIGURE 8

FIGURE 8. The 13 elastic stiffness constants from the monoclinic model plotted against the two sets of ${\mathit{\varepsilon }}_{\mathit{1}}$ , ${\mathit{\varepsilon }}_{\mathit{2}}$ while ${\psi }_{\mathit{1}},{\psi }_{\mathit{2}}$ for two sets are fixed at 35^° and 45^° respectively.

In Figure 9, the 13 elastic stiffness constants have been plotted as a function of ${\psi }_{\mathit{1}}$ and ${\psi }_{\mathit{2}}$ for the two fracture sets with constant ${\epsilon }_{\mathit{1}}$ and ${\epsilon }_{\mathit{2}}$ at 0.03 and 0.05, respectively. There is a periodic trend of each elastic stiffness constant with ${\psi }_{\mathit{1}}$ and ${\psi }_{\mathit{2}}$ , thus following the transformation law of the fourth ranked tensor given by (Auld, 1990). There will be no periodic trend if fractures are not present within the formation.

FIGURE 9

FIGURE 9. The 13 elastic stiffness constants obtained using monoclinic model plotted as a function of ${\psi }_{\mathit{1}},{\psi }_{\mathit{2}}$ for two fracture sets. ${\mathit{\varepsilon }}_{\mathit{1}}$ , ${\mathit{\varepsilon }}_{\mathit{2}}$ for each set are kept constant in this case at 0.03 and 0.05.

4 Discussion

4.1 Seismic modeling using rock physics model for monoclinic symmetry

Once the saturated effective elastic properties of the porous fracture media have been calculated, the multiple fracture sets can then be characterized on the basis of seismic modeling. The variation of fracture based elastic properties within a rock can be detected based on sophisticated seismic attributes like seismic AVAZ, azimuthal variation of velocity in fractured interval, and shear wave birefringence analysis (Crampin et al., 1980; Lynn et al., 1995; Lynn et al., 1999; Zhu et al., 2004; Will et al., 2005). In this study, a practical method is presented for the development of rock physics model in case of a media containing multiple fracture sets utilizing conventional wireline log and 3D seismic post-stack data. The output from this rock physics model in terms of 13 independent effective elastic stiffness constants can be utilized for determination of monoclinic reflectivity given by Schoenberg and Protazio (1992):

$\mathbf{R}=\left[\begin{array}{ccc}{R}_{PP}& R_{SP}& R_{TP}\\ R_{PS}& R_{SS}& R_{TS}\\ R_{PT}& R_{ST}& R_{TT}\end{array}\right](11)$

$\mathit{X}=\left[\begin{array}{ccc}{e}_{P\mathit{1}}& e_{S\mathit{1}}& e_{T\mathit{1}}\\ e_{P\mathit{2}}& e_{S\mathit{2}}& e_{T\mathit{2}}\\ \{-\left(C_{\mathit{13}}{e}_{P\mathit{1}}+C_{\mathit{36}}{e}_{P\mathit{2}}\right)s_{\mathit{1}}& \{-\left(C_{\mathit{13}}{e}_{S\mathit{1}}+C_{\mathit{36}}{e}_{S\mathit{2}}\right)s_{\mathit{1}}& \{-\left(C_{\mathit{13}}{e}_{T\mathit{1}}+C_{\mathit{36}}{e}_{T\mathit{2}}\right)s_{\mathit{1}}\\ -\left(C_{\mathit{23}}{e}_{P\mathit{2}}+C_{\mathit{36}}{e}_{P\mathit{1}}\right)s_{\mathit{2}}& -\left(C_{\mathit{23}}{e}_{S\mathit{2}}+C_{\mathit{36}}{e}_{S\mathit{1}}\right)s_{\mathit{2}}& -\left(C_{\mathit{23}}{e}_{T\mathit{2}}+C_{\mathit{36}}{e}_{T\mathit{1}}\right)s_{\mathit{2}}\\ -C_{\mathit{33}}{e}_{P\mathit{3}}{s}_{\mathit{3}P}\}& -C_{\mathit{33}}{e}_{S\mathit{3}}{s}_{\mathit{3}S}\}& -C_{\mathit{33}}{e}_{T\mathit{3}}{s}_{\mathit{3}T}\}\end{array}\right](12)$

$\mathit{Y}=\left[\begin{array}{ccc}\{-\left(C_{\mathit{55}}{s}_{\mathit{1}}+C_{\mathit{45}}{s}_{\mathit{2}}\right)e_{P\mathit{3}}& \{-\left(C_{\mathit{55}}{s}_{\mathit{1}}+C_{\mathit{45}}{s}_{\mathit{2}}\right)e_{S\mathit{3}}& \{-\left(C_{\mathit{55}}{s}_{\mathit{1}}+C_{\mathit{45}}{s}_{\mathit{2}}\right)e_{T\mathit{3}}\\ -\left(C_{\mathit{55}}{e}_{P\mathit{1}}+C_{\mathit{45}}{e}_{P\mathit{2}}\right)s_{\mathit{3}P}\}& -\left(C_{\mathit{55}}{e}_{S\mathit{1}}+C_{\mathit{45}}{e}_{S\mathit{2}}\right)s_{\mathit{3}S}\}& -\left(C_{\mathit{55}}{e}_{T\mathit{1}}+C_{\mathit{45}}{e}_{T\mathit{2}}\right)s_{\mathit{3}T}\}\\ & & \\ \{-\left(C_{\mathit{45}}{s}_{\mathit{1}}+C_{\mathit{44}}{s}_{\mathit{2}}\right)e_{P\mathit{3}}& \{-\left(C_{\mathit{45}}{s}_{\mathit{1}}+C_{\mathit{44}}{s}_{\mathit{2}}\right)e_{S\mathit{3}}& \{-\left(C_{\mathit{45}}{s}_{\mathit{1}}+C_{\mathit{44}}{s}_{\mathit{2}}\right)e_{T\mathit{3}}\\ -\left(C_{\mathit{45}}{e}_{P\mathit{1}}+C_{\mathit{44}}{e}_{P\mathit{2}}\right)s_{\mathit{3}P}\}& -\left(C_{\mathit{45}}{e}_{S\mathit{1}}+C_{\mathit{44}}{e}_{S\mathit{2}}\right)s_{\mathit{3}S}\}& -\left(C_{\mathit{45}}{e}_{T\mathit{1}}+C_{\mathit{44}}{e}_{T\mathit{2}}\right)s_{\mathit{3}T}\}\\ & & \\ e_{P\mathit{3}}& e_{S\mathit{1}}& e_{T\mathit{1}}\end{array}\right](13)$

Here, $\mathit{R}$ is reflection matrix; $s_{\mathit{1}}$ and $s_{\mathit{2}}$ defines the phase slowness vector (horizontal components), the associated eigenvectors are denoted by $e_P$ _, $e_S$_, and $e_T$ derived from the Christoffel equations (Mavko et al., 2009), and $C_{ij}$ denotes the elastic stiffness constants. For the case, when $\mathit{X}$ and $\mathit{Y}$ both are not singular and $\left({\mathit{X}}^{-\mathit{1}}\mathit{X}\mathit{ʹ}+{\mathit{Y}}^{-\mathit{1}}\mathit{Y}\mathit{ʹ}\right)$ can be mathematically inverted, the reflection matrix is given as:

$\mathit{R}=\left({\mathit{X}}^{-\mathit{1}}\mathit{X}\mathit{ʹ}-{\mathit{Y}}^{-\mathit{1}}\mathit{Y}\mathit{ʹ}\right){\left({\mathit{X}}^{-\mathit{1}}\mathit{X}\mathit{ʹ}+{\mathit{Y}}^{-\mathit{1}}\mathit{Y}\mathit{ʹ}\right)}^{-\mathit{1}}(14)$

Here, $\mathit{X}\mathit{ʹ}$ and $\mathit{Y}\mathit{ʹ}$ are similar to $\mathit{X}$ and $\mathit{Y}$ with only difference being replacement of unprimed parameters (incidence medium) by primed parameters (transmission medium). Once, reflectivity from above Eqs. (11-14) is obtained, it can be convolved with the seismic wavelet to obtain the synthetic (calculated) seismic AVAZ data at reservoir level.

4.2 Seismic inversion for porous media with multiple fracture sets

Inverse modeling targets the spatial estimation of fracture parameters, which influence the fluid flow during production (Ali et al., 2015). To obtain the modeled fracture parameters on a reservoir scale, the calculated seismic AVAZ data has to be matched with the field based observed seismic AVAZ data in order to minimize the error before the parameters under observation can be studied. The inverse problem in this particular case can be formulated as given in Eq. 15.

$G\left(\mathit{m}\right)\approx \mathit{d}.(15)$

Here, $\mathit{d}$ is a vector of observable quantities (acquired AVAZ data), $G$ combines the monoclinic rock physics model (as discussed in Rock Physics Modeling for Sakesar Formation) with the seismic model (as discussed in Seismic modeling using rock physics model for monoclinic symmetry), and $\mathit{m}$ inculcates the fracture parameters, selected for the developed monoclinic model (${\epsilon }_{\mathit{1}},{\epsilon }_{\mathit{2}}$ and ${\psi }_{\mathit{1}}$ , ${\psi }_{\mathit{2}}$), in the form of a vector. The workflow for inversion of parameters using correlation of synthetic AVAZ data with acquired AVAZ data is given in Figure 10.

FIGURE 10

FIGURE 10. Workflow proposed by the study to invert for fracture parameters of multiple fracture sets by developing synthetic AVAZ volume using rock physics model for monoclinic symmetry which can be correlated with observed 3D seismic AVAZ data.

Correlation of the developed synthetic data with the pre-stack 3D seismic AVAZ data will lead to spataial mapping the parameters at reservoir scale, thus immensely assisting in reservoir characterization provided that the synthetic data has been modelled correctly. The plots for these parameters with 13 effective elastic stiffness constants helps to monitor the sensitivity of the forward model. Changing the forward model will change the outcome of the reflectivity. Therefore, the geostatistical correlation helps to select the best possible model to invert for the fracture parameters.

Due to limited availability of borehole data over large areas, a monoclinic rock physics model can be opted to spatially characterize the azimuthal fracture orientations and fracture densities of multiple fracture sets in a fractured reservoir. The challenge however is requirement of accurate priori information regarding petrophysical properties like fluid saturation, porosity, volume of shale etc.; matrix properties like matrix moduli, matrix density, and fracture geometry. The importance of rock physics model lies in its ability to inculcate a large number of input parameters. This however can turn out to be a problem provided an improper priori information is used. It is therefore essential to know the alignment of fracture sets present in a fractured reservoir before a rock physics model can be developed. The Ant-tracking DFN modeling proves to be a good tool for this purpose, however, core data analysis is always preferable. The aspect ratio for fractures used in this study was obtained from outcrop of the studied formation, but core sample analysis for aspect ratio determination is suggested as it affects the results of T-matrix model.

5 Conclusion

In the absence of sophisticated laboratory/wireline data, characterization of fracture parameters (like azimuthal orientation and density) of multiple fracture sets in a fractured media can be achieved using a suitable rock physics model. The accuracy of the developed model depends upon the precision of input parameters and awareness of background geological knowledge. The Ant-tracking DFN model gives a potential initial estimate of fractures, which can be used as a priori information for the development of the required rock physics model. Here, in this study, a proper rock physics model for a media comprising of multiple fracture sets (monoclinic symmetry) was developed using T-matrix approach, Wood’s equation, and Brown-Korringa relation. The output in the form of 13 independent effective elastic stiffness coefficients shows an inverse linear trend with ${\epsilon }_{\mathit{1}}$ and ${\epsilon }_{\mathit{2}}$ There is a periodic variation of the effective stiffness coefficients with ${\psi }_{\mathit{1}}$ and ${\psi }_{\mathit{2}}$ confirming the presence of fractures within the formation. These elastic constants can practically be applied to generate a synthetic pre-stack (calculated) seismic AVAZ data using the Schoenberg and Protazio (1992) solution. Fracture properties can effectively be inverted through geostatistical correlation of this synthetic data with pre-stack (observed) 3D seismic AVAZ data. This research facilitates and provides a practical approach to develop a proper rock physics model applicable to fractured reservoirs in order to study the parameters of existing fractures, thus facilitating in increasing the production from these reservoirs.

Author's note

This work is a part of FM's Ph.D thesis.

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: The data was officially obtained after applying through the university to the Directorate General of Petroleum Concessions (DGPC) which has jurisdiction over handling the data. Requests to access these datasets should be directed to dgpc@mpnr.gov.pk.

Author contributions

FM implemented the methodology and carried out the research. ZA proof read the manuscript and contributed in finalizing the manuscript. YA provided technical support and assisted in carrying out the methodology. Zubair assisted in specialized technical support. AA put forward the idea and helped with the provision of methodology for this research as well as assessed the results of the study.

Acknowledgments

We would like to thank the Directorate General of Petroleum Concessions (DGPC), LMKR, Pakistan, for being a data source for this work and the Department of Earth Sciences, Quaid-i-Azam University (QAU), Islamabad, Pakistan, which provided the basic facilities for completing this work. We personally thank Tahir Hussain, Ph.D. Scholar, QAU, Muhammad Zahid, Chief Petrophysicist, Mari Petroleum Company Limited, and Ali Zulqarnain, Geophysicist, OGDCL, who provided help and support during this case study.

Conflict of interest

Author ZA is employed by the Digital and Integration, Schlumberger, Pakistan

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Ali, A., and Jakobsen, M. (2011a). On the accuracy of rüger’s approximation for reflection coefficients in HTI media: Implications for the determination of fracture density and orientation from seismic AVAZ data. J. Geophys. Eng. 8, 372–393. doi:10.1088/1742-2132/8/2/022

CrossRef Full Text | Google Scholar

Ali, A., Anwer, H. M., and Hussain, M. (2015). A comparisonal study in the context of seismic fracture characterization based on effective stiffness and compliance methods. Arab. J. Geosci. 8, 4117–4125. doi:10.1007/s12517-014-1478-8

CrossRef Full Text | Google Scholar

Ali, A., and Jakobsen, M. (2014). Anisotropic permeability in fractured reservoirs from frequency-dependent seismic Amplitude versus Angle and Azimuth data. Geophys. Prospect. 62 (2), 293–314. doi:10.1111/1365-2478.12084

CrossRef Full Text | Google Scholar

Ali, A., and Jakobsen, M. (2011b). Seismic characterization of reservoirs with multiple fracture sets using velocity and attenuation anisotropy data. J. Appl. Geophys. 75 (3), 590–602. doi:10.1016/j.jappgeo.2011.09.003

CrossRef Full Text | Google Scholar

Auld, B. A. (1990). Acoustic fields and waves in solids. Malabar, FL: Krieger Publishing.

Google Scholar

Aydin, H., and Akin, S. (2019). “Discrete fracture network modeling of alaşehir geothermal field,” in Proceedings of the 44th Workshop on Geothermal Reservoir Engineering, Stanford, California, 11-13 February (Stanford University). SGP-TR-214.

Google Scholar

Bagni, F. L., Bezerra, F. H., Balsamo, F., Maia, R. P., and Dall'Aglio, M. (2020). Karst dissolution along fracture corridors in an anticline hinge, jandaíra formation, Brazil: Implications for reservoir quality. Mar. Petroleum Geol. 115, 104249. doi:10.1016/j.marpetgeo.2020.104249

CrossRef Full Text | Google Scholar

Batzle, M., and Wang, Z. (1992). Seismic properties of pore fluids. Geophysics 57 (11), 1396–1408. doi:10.1190/1.1443207

CrossRef Full Text | Google Scholar

Bender, F., Raza, H. A., and Bannert, D. N. (1995). Geology of Pakistan. Berlin, Stuttgart: Gebruder Borntraeger.

Google Scholar

Bisdom, K., Gauthier, B. D. M., Bertotti, G., and Hardebol, N. J. (2014). Calibrating discrete fracture-network models with a carbonate three-dimensional outcrop fracture network: Implications for naturally fractured reservoir modeling. Am. Assoc. Pet. Geol. Bull. 98 (7), 1351–1376. doi:10.1306/02031413060

CrossRef Full Text | Google Scholar

Brown, R. J. S., and Korringa, J. (1975). On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics 40, 608–616. doi:10.1190/1.1440551

CrossRef Full Text | Google Scholar

Carcione, J. M. (1995). Constitutive model and wave equations for linear, viscoelastic, anisotropic media. Geophysics 60, 537–548. doi:10.1190/1.1443791

CrossRef Full Text | Google Scholar

Cox, T., and Seitz, K. (2007). “May. Ant tracking seismic volumes for automated fault interpretation,” in CSPG CSEG convention (Oklahoma, United States: AAPG), 670571.

Google Scholar

Crampin, S., McGonigle, R., and Bamford, D. (1980). Estimating crack parameters from observations of P-wave velocity anisotropy. Geophysics 45, 345–360. doi:10.1190/1.1441086

CrossRef Full Text | Google Scholar

Fang, J., Zhou, F., and Tang, Z. (2017). Discrete fracture network modelling in a naturally fractured carbonate reservoir in the jingbei oilfield, China. Energies 10 (2), 183. doi:10.3390/en10020183

CrossRef Full Text | Google Scholar

Gee, E. R., and Gee, D. G. (1989). Overview of the geology and structure of the Salt Range, with observations on related areas of northern Pakistan, 232. Boulder, Colorado: Geological Society of America Special Paper, 95–112.

CrossRef Full Text | Google Scholar

Gholipour, A. M., Cosgrove, J. W., and Ala, M. (2016). New theoretical model for predicting and modelling fractures in folded fractured reservoirs. Pet. Geosci. 22, 257–280. doi:10.1144/petgeo2013-055

CrossRef Full Text | Google Scholar

Hu, J. L., Kang, Z. H., and Yuan, L. L. (2014). “Automatic fracture identification using ant tracking in Tahe oilfield,” in Advanced materials research. Editors X. Lu, K-C. Lam, H. Xu, and Z. Jia (Freienbach, Switzerland: Trans Tech Publications Ltd), 962, 556–559.

CrossRef Full Text | Google Scholar

Huang, S., Zhang, Y., Zheng, X., Zhu, Q., Shao, G., Cao, Y., et al. (2017). Types and characteristics of carbonate reservoirs and their implication on hydrocarbon exploration: A case study from the eastern tarim basin, NW China. J. Nat. Gas Geoscience 2 (1), 73–79. doi:10.1016/j.jnggs.2017.02.001

CrossRef Full Text | Google Scholar

Jadoon, M. S. K., Jadoon, I. A. K., Bhatti, A. H., and Ali, A. (2006). Fracture characterization and their impact on the field development. Pak. J. Hydrocarbon Res. 16, 11–21. doi:10.2118/111053-MS

CrossRef Full Text | Google Scholar

Jakobsen, M., Hudson, J. A., and Johansen, T. A. (2003a). T-Matrix approach to shale acoustics. Geophys. J. Int. 154, 533–558. doi:10.1046/j.1365-246x.2003.01977.x

CrossRef Full Text | Google Scholar

Jakobsen, M., Johansen, T. A., and McCann, C. (2003b). The acoustic signature of fluid flow in complex porous media. J. Appl. Geophy. 54, 219–246. doi:10.1016/j.jappgeo.2002.11.004

CrossRef Full Text | Google Scholar

Kemal, A., Balkwill, H. R., Stoakes, F. A., Ahmad, G., Zaman, A. S. H., and Humayon, M. (1991). “Indus basin hydrocarbon plays,” in International Petroleum Seminar on new directions and strategies for accelerating Petroleum Exploration and Production in Pakistan, 16–57.

Google Scholar

Khan, M. A., Ahmed, R., Raza, H. A., and Kemal, A. (1986). Geology of petroleum in Kohat-Potwar depression, Pakistan. AAPG Bull. 70 (4), 396–414. doi:10.1306/9488571E-1704-11D7-8645000102C1865D

CrossRef Full Text | Google Scholar

Kumar, D. (2006). A tutorial on Gassmann fluid substitution: Formulation, algorithm and matlab code. Geohorizons 11, 4–12.

Google Scholar

Lynn, H. B., Beckham, W. E., Simon, K. M., Bates, C. R., Layman, M., and Jones, M. (1999). P-wave and S-wave azimuthal anisotropy at a naturally fractured gas reservoir, Bluebell-Altamont field, Utah. Geophysics 64, 1312–1328. doi:10.1190/1.1444636

CrossRef Full Text | Google Scholar

Lynn, H. B., Simon, K. M., Bates, C. R., Layman, M., Schneider, R., and Jones, M. (1995). Use of anisotropy in P-wave and S-wave data for fracture characterization in a naturally fractured gas reservoir. Lead. Edge 14, 887–893. doi:10.1190/1.1437179

CrossRef Full Text | Google Scholar

Masood, F., Ahmad, Z., and Khan, M. S. (2017). Moderate interpretation with attribute analysis and 3D visualization for deeper prospects of balkassar field, central potwar, upper Indus Basin, Pakistan. Int. J. Geosciences 8 (05), 678–692. doi:10.4236/ijg.2017.85037

CrossRef Full Text | Google Scholar

Mavko, G., Mukerji, T., and Dvorkin, J. (2009). The rock physics handbook. Cambridge, United Kingdom: Cambridge University Press.

Google Scholar

Misaghi, A., Negahban, S., Landrø, M., and Javaherian, A. (2010). A comparison of rock physics models for fluid substitution in carbonate rocks. Explor. Geophys. 41, 146–154. doi:10.1071/eg09035

CrossRef Full Text | Google Scholar

Narr, W., Schechter, D. S., and Thompson, L. B. (2006). Naturally fractured reservoir characterization, 112. Richardson, Texas: Society of Petroleum Engineers.

Google Scholar

Nguyen, P. K. T., and Nam, M. J. (2011). A review on methods for constructing rock physics model of saturated reservoir rock for time-lapse seismic. Geosystem Eng. 14 (2), 95–107. doi:10.1080/12269328.2011.10541336

CrossRef Full Text | Google Scholar

Ouenes, A. (2000). Practical application of fuzzy logic and neural networks to fractured reservoir characterization. Comput. Geosci. 26 (8), 953–962. doi:10.1016/s0098-3004(00)00031-5

CrossRef Full Text | Google Scholar

Pedersen, S. I., Randen, T., Sonneland, and L., and Steen, Ø. (2002). “Automatic fault extraction using artificial ants,” in SEG technical program expanded abstracts 2002 (Texas, United States: Society of Exploration Geophysicists), 512–515.

CrossRef Full Text | Google Scholar

Sayers, C. M. (2009). Seismic characterization of reservoirs containing multiple fracture sets. Geophys. Prospect. 57 (2), 187–192. doi:10.1111/j.1365-2478.2008.00766.x

CrossRef Full Text | Google Scholar

Schoenberg, M., and Protazio, J. (1992). Zoeppritz rationalized and generalized to anisotropy. J. Seism. Explor. 1, 125–144. doi:10.1121/1.2029011

CrossRef Full Text | Google Scholar

Silva, C. C., Marcolino, C. S., and Lima, F. D. (2005). “Automatic fault extraction using ant-tracking algorithm in the marlim south field, campos basin,” in SEG technical program expanded abstracts (Texas, United States: Society of Exploration Geophysicists), 857–860.

CrossRef Full Text | Google Scholar

Smith, T. M., Sondergeld, C. H., and Rai, C. S. (2003). Gassmann fluid substitutions: A tutorial. GEOPHYSICS 68 (2), 430–440. doi:10.1190/1.1567211

CrossRef Full Text | Google Scholar

Souque, C., Knipe, R. J., Davies, R. K., Jones, P., Welch, M. J., and Lorenz, J. (2019). Fracture corridors and fault reactivation: Example from the chalk, isle of thanet, kent, england. J. Struct. Geol. 122, 11–26. doi:10.1016/j.jsg.2018.12.004

CrossRef Full Text | Google Scholar

Tran, N. H., Chen, Z., and Rehman, S. S. (2006). Integrated conditional global optimisation for discrete fracture network modelling. Comput. Geosciences 32, 17–27. doi:10.1016/j.cageo.2005.03.019

CrossRef Full Text | Google Scholar

Wang, B., Chen, Y., Lu, J., and Jin, W. (2018). A rock physics modelling algorithm for simulating the elastic parameters of shale using well logging data. Sci. Rep. 8 (1), 12151–12158. doi:10.1038/s41598-018-29755-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Will, R., Archer, R., and Dershowitz, B. (2005). Integration of seismic anisotropy and reservoir performance data for characterization of naturally fractured reservoirs using discrete fracture network models, SPE Reserv. Eval. Eng. 8 (2), 132–142. doi:10.2118/84412-MS

CrossRef Full Text | Google Scholar

Wood, A. W. (1955). A textbook of sound. New York: Macmillan, 360.

Google Scholar

Zhiqian, G., Zhongbao, L., Shanlin, G., Qunan, D., Shiqiang, W., and Shilin, L. (2016). Characteristics and genetic models of Lower Ordovician carbonate reservoirs in southwest Tarim Basin, NW China. J. Petroleum Sci. Eng. 144, 99–112. doi:10.1016/j.petrol.2016.03.007

CrossRef Full Text | Google Scholar

Zhu, P., Wang, W., Yu, W., and Zhu, G. (2004). Inverting reservoir crack density from P-wave AVOA data. J. Geophys. Eng. 1, 168–175. doi:10.1088/1742-2132/1/2/010

CrossRef Full Text | Google Scholar

Zhu, X., and McMechan, G. A. (1990). Direct estimation of the bulk modulus of the frame in fluid saturated elastic medium by Biot theory. Seg. Tech. Program Expand. Abstr. 1990, 787–790. doi:10.1190/1.1890340

CrossRef Full Text | Google Scholar