1. Introduction

There are several techniques for estimating the elastic properties of the earth’s subsurface; these are also used for other purposes such as reservoir characterization by means of rock physics petro-elastic models (PEM) or/and fluid and lithology detection. To this end, amplitude-variation-with-offset (AVO) analysis and inversion of the seismic data have been frequently used by geophysicists for many years. In this connection, AVO inversion is an area of research that attempts to go beyond the qualitative process of AVO interpretation used to predict the presence of gas in a sand (Hampson, 1991). The following is a brief explanation of the inversion methods applied in this context.

Post-stack inversion methods are the first approaches for inverting the seismic data into elastics properties. The early methods of post-stack inversion were recursive methods that suffered from recovering the low-frequency component of the reflectivity. After that, many inversion methods such as colored inversion method with spectral bluing (Neep, 2007) and model-based inversion were developed to solve the problems of the recursive method (Cooke and Cant, 2010).

One of these methods for estimating the elastic properties of the earth is elastic impedance, which was introduced by Connolly (1999) and improved by Whitcombe et al. (2002) as an extended elastic impedance (EEI) generalized form of the acoustic impedance. This method first inverts angle-limited stacks under the constraint of simulated logs and then extracts the elastic properties. After extracting the elastic properties, these data can be converted to the interesting information such as delineating hydrocarbon bearing zones within a reservoir (Karbalaali et al., 2013), porosity, shale content, and water saturation (Samba et al., 2017). These types of inversion (post-stack and EEI inversion) and reservoir characterization are not studied in this paper.

Alternatively, AVO inversion is another method of inversion that includes extracting AVO attributes from seismic data and then inverting them under the constraint of the appropriate logs (Cambois, 2000). The amplitude reflection of P-wave can be approximated by many relations that contain AVO attributes such as those introduced by Shuey (1985), Fatti et al. (1994), and Verm and Hilterman (1995) that include zero offset P-wave reflectivity (R_P), zero offset S-wave reflectivity (R_S), AVO intercept (A), AVO gradient (B), and Poisson’s reflectivity (PR). In addition to these attributes, Goodway et al. (1997) and Veeken and Rauch-Davies (2006) introduced more AVO attributes which are known as Lambda-Mu-Rho (LMR). In fact, these attributes are the results of multiplication of the Lame parameters and density. Then, under the constraint of logs, these attributes were converted to obtain the acoustic impedance (P and S-wave), V_P/V_S, Lambda-Rho, and Mu-Rho. Because AVO attributes are sensitive to the presence of noise such as wavelet variation with offset, poor normal moveout (NMO) corrections, and residual noise (multiples and converted waves), this method of inversion may provide unreliable results. In addition to the above method, there are other ways of AVO inversion (pre-stack inversions, in general) such as those introduced by Simmons and Backus (1996) and Buland and Omre (2003). In the first method, Simmons and Backus (1996), inversion is done for linearized P-reflectivity, S-reflectivity, and density regarding the following three assumptions:

• The reflectivity terms can be estimated from the angle-dependent reflectivity by Aki and Richards linearized approximation.
• Density and P-wave velocity are related by Gardner’s equation (Gardner et al., 1985).
• S-wave velocity and P-wave velocity are related by mudrock line (Castagna et al., 1985).

Therefore, the reflectivity is given using the linearized inversion approach. The second method uses a similar inversion technique, i.e. Bayesian linearized AVO inversion. Their method is parameterized by three terms: P-reflectivity, S-reflectivity, and ρ reflectivity. Then, using the Aki and Richards’ approximation and the assumption of the small reflectivity approximation, the changes of these parameters are related to the original parameter itself. The small reflectivity assumption allows the user to perform the inversion for velocity and density, rather than reflectivity. Hampson et al. (2005) developed a new approach, namely simultaneous inversion, that allows the user to invert directly P-impedance, S-impedance, and density. Avseth et al. (2016) implemented this method in the Norwegian Sea for the hydrocarbon exploration; they pointed out that considering geology information of the area is necessary for such studies. This type of inversion is considered as a deterministic approach. There are two major statistical frameworks for AVO inversion; i.e. frequentist framework and Bayesian framework. The first method assumes that there is a fixed but unknown truth about the subsurface and seeks to estimate that truth as closely as possible. By contrast, the second method describes the understanding of the subsurface as a prior probability distribution. The inversion approach does not involve the modification of the forward model to make it more stable. Instead, those data are used in combination with the forward model to produce a posterior probability distribution that characterizes understanding of the subsurface given the seismic data (Ball et al., 2016).

In the present study, simultaneous inversion method is performed on one of the oilfields located on the western of the Persian Gulf. The reservoir target is sandstone with the interlayers of shale and dolomite.

2. Simultaneous inversion of pre-stack seismic data

Inversion is a method for inverting amplitude (one of the components of the seismic data) of the seismic reflection data for many purposes. The other component of the seismic reflection data is the travel time, which is used for constructing the geometry of the subsurface for hydrocarbon industry. The first one is known as stratigraphic inversion, and the second one is known as structure inversion (Yilmaz, 2008). The inversion methods are divided into two groups, which use post-stack seismic data and pre-stack data. The post-stack techniques transform the data to density and compressional wave velocity, which is known as P-impedance among the geophysicists. Then, using the inverted data, geophysicists are able to make predictions about some petrophysical and reservoir properties such as porosity and lithology. However, these predictions are somehow ambiguous because P-impedance is affected by many parameters such as lithology, porosity, and pore fluid; moreover, it is difficult to separate their effects from each other. Although post-stack inversion methods are quick, they lack the separation of the effects of the porosity, pore fluids, and lithology. In order to make a comparison, a model-based inversion and simultaneous inversion result are shown in Figure 1, which shows the limitation of model-based inversion for lithology discrimination. ‌The capability of this method will be discussed in the following sections. To improve these ambiguous results, the full elastic inversion, in which P-impedance, S-impedance, and density are estimated, is performed (Russell and Hampson, 2006). This type of inversion falls in the second group, which is known as pre-stack inversion method.

Figure 1

a) Model-based inversion result and b) simultaneous inversion. The color bars show P-impedance values ((m/s)×(g/cc)); the vertical black line is the well path and the black log is V_p. Section a) is brought from crossline 1812 and inline 146 to 284 and section b) is brought from crossline 1812 and inline 230 to 272.

The goal and the purpose of the pre-stack seismic inversion are to estimates P-impedance, S-impedance, and density, followed by predicting reservoir properties such as the property of fluid and rock. In other words, pre-stack inversion can be considered as a quantitative extension of AVO (Russell and Hampson, 2006). Hence, the approximated relations of Zoeppritz are needed. The re-formulated Aki and Richards’ equation is given below (Fatti et al., 1994)

where, R_P is zero offset P-wave reflectivity, and R_S represents zero offset S-wave reflectivity; R_D is also density reflectivity.

Furthermore, the linearized zero offset P-reflectivity, S-reflectivity, and density reflectivity are given by (Aki and Richards, 2002):

where, , , and  are difference of P-wave, S-wave velocity, and density of the two-layer structure respectively, and V_p, V_s, and ρ are the averages of P-wave, S-wave velocity, and density of the two-layer structure respectively. After extracting the reflectivity in Equation 1, they can be converted by post-stack inversion methods such as model-based inversion. This is referred to as independent inversion (Russell and Hampson, 2006). Hampson et al. (2005) developed a new algorithm by combining algorithm of Simmons and Backus (1993) and Buland and Omre (2003). This algorithm is based on three assumptions (Hampson et al., 2005):

• The linearized approximation for reflectivity holds;
• PP and PS reflectivity as a function of angle can be given by the Aki and Richards equations;
• The relationship between the logarithm of P-impedance and both S-impedance and density is linear.

Therefore, their algorithm uses angle gather data and allows inverting directly for P-impedance, S-impedance, and density. This method is referred to as simultaneous inversion. They first developed it for post-stack data and then extended it for pre-stack seismic data. The small reflectivity approximation of the P-wave reflectivity is given by (Hampson et al., 2005)

where, i represents the i^th interface between layers i and i+1, and Z_P is P-impedance. Thus, Equation 5 can be written in a matrix form as follows (Hampson et al., 2005):

where,  Then, if the seismic trace can be considered as a convolution of the seismic wavelet, the Earth’s reflectivity is expressed by:

where, T_i represents the i^th sample of the seismic trace, and W_j is the j^th term of the extracted seismic wavelet; * stands for convolution. Combining Equations 6 and 7 gives a relationship between seismic trace and logarithm of P-impedance as defined by (Hampson et al., 2005):

where, W is wavelet matrix, and D is derivative matrix given in Equation 6. By combining Equations 1 and 8 for pre-stack case (Hampson et al., 2005), we have:

where,  and  Moreover, there is a linear relation between L_S, L_P, and L_D as follows (Hampson et al., 2005)

By combining Equations 9, 10, and 11, we obtain the following equation (Hampson et al., 2005):

where,  and The matrix form of Equation 12 is given by:

It is noteworthy that  and  are a deviation from a straight line when ln(Z_S) is plotted versus ln(Z_P); ln(ρ) is plotted versus ln(Z_P) respectively. The deviations from the straight line,  and , are the desired fluid anomalies (Hampson et al., 2005); eventually, we solve Equation 13. Using the standard matrix inversion technique for this purpose, two problems arise. The first one is the cost and the unstable potential of the matrix inversion, and the second one, which is more important, is the low-frequency component that will not recover. Another method for solving this problem is to build an initial P-impedance model, and then iterate toward a solution using the method of conjugate gradients, which is recommended. Finally, the results of inversion are obtained.

3. AVO inversion workflow

This method of AVO inversion was performed on one of the western oilfields of the Persian Gulf, which is located in the south of Iran. The reservoir consists of Oligo-Miocene unconsolidated sands, with interbedded shale, dolomite, dolomite cemented sandstones, and layers of nodular anhydrite, and the base of the reservoir is composed of carbonate rocks. The field has almost a constant thickness and is capped by a few meters of anhydrite, which makes the base of the Upper Asmari formation. The reservoir at the crest of the structure is at a depth of 815 m with a hydrocarbon column of about 62 m, of which the gas column is 18 m, and oil column is 44 m with an oil gravity of 27 API (Ghazban, 2009). The sands of the Ghar member are fine- to medium-grained with some shale and dolomite interbeds. The Ghar member has excellent reservoir properties with porosity in the order of 35% and permeability of around 2 to 4 mD. In addition, the sands are clean and poorly cemented. The field is a NW-SE trending anticline, dipping by 1-1.5 degrees measuring at approximately 24 by 8 km, and filled to spill point by a hydrocarbon column of about 62 m (Ghazban, 2009). The structure is an elongated and very gently deformed NW-SE symmetrical anticline which is parallel to the direction of the Zagros mountains.

The methodology of this study is summarized in Figure 2. There are 3 wells in the studied area which have P-wave velocity, density, and gamma ray. Only one of them has S-wave velocity measurement, so the estimate of S-wave velocity is required for the rest of them. Therefore, the S-wave velocity was calculated by means of Xu and White’s (1995) relation. In fact, the estimates of S-wave velocity may be obtained from known mineral matrix properties and measured porosity and clay content or from measured P-wave velocity and either porosity or clay content (Mavko et al., 2009). An example of the logs is illustrated in Figure 3. Moreover, check shot corrections were used to correct the P-wave velocity logs. The processing steps of seismic pre-stack data used for inversion purposes were statics corrections, noise attenuation, deconvolution, demultiplying, amplitude recovery, NMO corrections, and velocity analysis. Moreover, we applied the amplitude recovery step for recovering the true amplitudes such as geometrical spreading and adsorption. Table 1 shows the processing scheme used for seismic data of this field. Furthermore, Ramos et al. (1999) proposed a method that was designed to properly deal with geometrical spreading losses with offset. This method is known as true amplitude DMO. After the application of this method, geometrical spreading losses are automatically transformed and reflection coefficients remain unchanged. Furthermore, this method is useful for area with small structural complexity such as the investigated field in the current paper. The main advantage of true amplitude DMO, in addition to reducing geometrical spreading losses, is a better estimation of AVO attributes such as intercept and gradient. Figure 4 shows a comparison between conventional AVO gradient analysis and true amplitude DMO. A considerably smaller scatter of the amplitudes around the background trend is observed due to poor gradient estimation. However, the cross plot of the true amplitude DMO shows a significantly larger scatter and better separation due to better geometrical spreading compensation of amplitudes at far offsets.

Figure 2

The workflow of the simultaneous AVO inversion.

Figure 3

Well data from one of the wells for reservoir section; from left, the curves are gamma ray, P-wave, S-wave, density, total porosity, and shale volume. 

Figure 4

Comparison between a) conventional AVO gradient analysis and b) true amplitude DMO.

Table 1

The processing scheme of the seismic data.

Seismic data quality plays an important role in data analysis and directly effects the results of inversion process. The data used in this study is 40-fold pre-stack seismic data with the dominant frequency of 35 Hz. Seismic data quality of the reservoir is superior because the reservoir is not deep, and the amplitude would be relatively preserved as opposed to deep reservoirs. For AVO inversion, amplitude recovery steps are very important and should be conducted carefully to show the true changes of amplitude with offset or angle. On this field, Ghar member is at a shallow depth and sufficiently thick to be resolved seismically, and the seismic data is of excellent quality. Therefore, it is an excellent case for seismic inversion studies. The Ghar member sandstone reservoir is characterized by low acoustic impedance through the entire section. Additionally, the gas production zone is characterized by low acoustic impedance. The dolomite or dolomite cemented sand interbeds are defined by higher acoustic impedance.

After the quality-controlling step, a velocity model was prepared for converting offset domain into an angle domain (i.e. angle gather). Since the simultaneous inversion algorithm needed wavelet in the angle domain, two statistical wavelets from two angle ranges were extracted from the resulted angle gathers. As expected, by increasing the angle, the wavelet was attenuated (Figure 5). Then, by means of the extracted wavelets, well tie analysis was performed by constructing synthetic data. In fact, this process was carried out by generating synthetic seismic data near the well location and correlating the synthetic data with the field seismic data. This kind of well to seismic tie method is common in the inversion analysis, so Figure 5c shows this process for one of the wells in the studied area. As shown in Figure 6, inversion analysis was applied to finding the constant parameters in Equations 10 and 11. The final step of this methodology is iterating Equation 13 toward a solution using the method of conjugate gradients and obtaining the inversion results. Figure 7 shows inversion results such as S-impedance and density. Alternatively, the velocities of the rock for a more accurate analysis were inverted (Figure 8a). Therefore, it is suitable to plot the cross-plot of P-impedance with respect to V_P/V_Sto determine lithology and fluid indicators. By means of these inverted elastic properties, one can also convert them to Lambda-Mu-Rho (LMR) parameters for better fluid detection. In the LMR analysis, the Lambda-Rho is a good attribute for this purpose. This fluid detection for sandstone member is illustrated in Figure 8b.

Figure 5

Statistical extracted wavelets from the near angle (red) and from the far angle (green); a) in the time domain and b) in the frequency domain; as expected, wavelets are changed by increasing angle (offset); c) well tie analysis was done by constructing synthetic data.

Figure 6

Inversion analysis for estimating the constants in Equations 10 and 11; a) cross-plot of ZS versus ZP and b) cross-plot of density versus ZP. The deviations away from the straight line,  and , are the desired fluid anomalies. The color bar shows the true vertical depth from SRD.

Figure 7

Inversion results; a) inverted S-impedance for the desired reservoir and b) inverted density; the vertical black line is the well path, and the black ellipse is dolomite cemented sandstone; the sections are from crossline 1812.

Figure 8

Fluids contact detection; a) V_p/V_s ratio and b) Lambda-Rho; the black curve XX’ is the WOC contact, and the black vertical line is the well path; the sections are from crossline 1812.

4. Results and discussion

The main purpose of AVO inversion is to estimate the subsurface rock properties using seismic P-P reflection pre-stack data as the input. As mentioned before, the reservoir studies in this work, which is located between top Ghar and top lower Asmari, is a sandstone gas and oil reservoir that includes interlayers of shale and dolomite. The presence of gas and oil results in some changes in hydrocarbon fluid and elasticity properties of the reservoir and consequently affects the seismic response of the reservoir. In addition to the effect of fluid on the elastic properties of rock, the minerals of rock, porosity, framework, and lithology could also influence the elasticity properties and therefore the seismic response. Accordingly, by the occurrence of these events, seismic can detect these changes. Moreover, gas oil contact (GOC) and oil-water contact (OWC) have similar effects which could be detected. As presented in Figure 7b, the reservoir zone is divided into five sections, which allow detecting the interlayers. It should be noted that this inversion is sensitive to lithology changes. As shown in Figure 7a, the presence of hydrocarbon (especially gas) decreases the inverted P-impedance. From the viewpoint of detection ability, the P-impedance, on its own, is not a gas sand or in general hydrocarbon indicator. However, by combining the ratio of the inverted P- and S-impedances, which gives the V_P/V_S ratio, the interpretation becomes clearer and more efficient. Along with a drop in P-impedance, the V_P/V_S ratio also declines, which is generally an indicator of a hydrocarbon and lithology. As expected, the V_P/V_S ratio is more sensitive and suitable for the purpose of differentiation of the fluid and lithology (Figure 8a). Furthermore, S-wave velocity has the main role in this differentiation because of its nature. Since the density of the hydrocarbon and water are different, the WOC could be detected for one of the wells in the studied area as displayed in Figure 8a, and by means of Lambda-Rho, this fluid contact is obvious (Figure 8b). As expected, the differences between V_P/V_Sand Lambda-Rho is a key for fluid contact detection. A dolomite-cemented sandstone interlayer is also seen in Figure 7b (black ellipse).

The geological surveys of the Ghar reservoir confirm that the reservoir thickness is constant, and it contains unconsolidated sand with interlayered shale, dolomite, dolomite cemented sandstone, and layers of nodular anhydrite. The cap rock of this reservoir is the thin anhydrite layer of the Upper Asmari formation; the base of the Ghar reservoir consists of carbonates. The results of AVO inversion are consistent with the geology surveys of the area. The low V_P/V_S ratio of the gas cap and the values of this ratio are reasonable for the base of the Ghar formation, which consists of carbonates. To show this result more conclusively, the cross-plot of the P-impedance and V_P/V_S ratio is required. At first, we examine this procedure on a well and check it with the geological reports, and we then use it for interpreting the inversion results. Therefore, Figure 9 shows the first step mentioned above, and Figure 10 displays the interpretation of the inversion results. The capability of velocity ratio is clear since the hydrocarbon-bearing sand has a lower P-impedance and V_P/V_S ratio, and the above carbonate formation has higher values (Figure 10). The main result of this inversion is providing the cubes of elastic properties of the formation or the reservoir rock such as V_P/V_S ratio, which then, by applying a good rock physics model, leads to a clear and suitable description of the reservoir (i.e. reservoir characterization). Once V_P, V_S, and velocities ratio cubes are calculated from the inversion results, by some additional calculations, rock elastic properties can be calculated. According to the relations between the rock properties and the rock elastic properties in an isotropic and homogenous medium, two sections of Poisson’s ratio and Young’s modulus are calculated (Figure 11). However, for a more accurate analysis, it is recommended that stochastic inversion methods be applied.

Figure 9

a) cross-plot of velocities ratio versus P-impedance for one of the wells in the studied area in which the red, yellow, and green colors show hydrocarbon-bearing sand, sandstone, and dolomite respectively; b) selected zones from a) are shown which are consistent with geological report.

Figure 10

a) cross-plot of V_P/V_S versus P-impedance, in which the red, yellow, and green zones show the hydrocarbon zones, sandstone, and dolomite zones respectively; b) the areas that are shown by the shapes (compare with Figure 9).

Figure 11

a) calculated Poisson’s ratio by the inversion results and b) calculated Young’s modulus (GPa) by the inversion results; both sections are from cross-line 1812 and the vertical black line shows the well path.

According to the results of the current study, AVO inversion has an excellent potential to reveal hydrocarbon reservoirs. However, there is a wide range of uncertainties about AVO inversion which can be grouped into three different categories: geological uncertainties, data uncertainties, and inversion uncertainties; hence, these uncertainties can be considered as the limitations of AVO inversion. The geological uncertainties affecting the inversion include compaction, overpressure, reservoir heterogeneity, fractures, and unknown geologies. There are uncertainties related to the acquisition, processing, and wave-propagation effects known as data uncertainties. Additionally, several uncertainties are associated with the inversion procedure such as the non-uniqueness of the inversion, which means there are several solutions which can provide the same seismic response; the need for a prior low-frequency starting model also represents an uncertainty during this kind of inversion. Another limitation to this kind of inversion is the assumption of an elastic and isotropic medium and the use of the three terms of Aki and Richards (2002) approximation to the Zoeppritz equations during AVO inversion.

5. Conclusions

In this paper, methods of estimation of elastic properties of subsurface formation are briefly discussed. However, only the simultaneous inversion is performed on the real data, and the results are examined. Our results show that P-impedance, on its own, is not sufficient and robust to discriminate between lithology and fluid; thus, the simultaneous inversion method was introduced. The AVO inversion of pre-stack seismic data by this method, constrained by geologic information, can greatly reduce interpretation risk. Moreover, this method of inversion provides an excellent parameter for discrimination study, i.e. V_P/V_S. Another main and robust use of AVO inversion is preparing a dataset for reservoir characterization, which is vital for the reservoir geophysics or generally reservoir management. By means of a good rock physics model and combining it with the inversion results, a good reservoir characterization is obtained, i.e. the cube of porosity. Although simultaneous AVO inversion is a great tool having the potential to reveal commercial hydrocarbon reservoir descriptions from pre-stack seismic data, there are several pitfalls and often large uncertainties associated with AVO inversion due to the deterministic nature of this method. However, this type of AVO inversion agrees with geological observations.

Nomenclature