Document Type: Research Paper
Authors
^{1} M.S. Student, Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Iran
^{2} 1Professor, Department of Petroleum Engineering, Amirkabir University of Technology, Tehran, Iran 2Professor, Institute of Geophysics, University of Tehran, Tehran, Iran
Abstract
Keywords
Main Subjects
Abstract Reservoir characterization has a leading role in the reservoir geophysics and reservoir management. Since the interests of the reservoir geophysics and reservoir managements are the elastic properties and reservoir properties of the subsurface rock for their purposes, a robust method is required for converting seismic data into elastic properties. Accordingly, by employing a rock physics model and using the inverted seismic data, one can describe the reservoir for purposes such as improvement in the production of the reservoir. In the present study, we employ one of the methods for converting the seismic data into the elastic properties. This method of inversion is known as simultaneous inversion, which is grouped in amplitudevariationwithoffset (AVO) inversion category. In this method, unlike the other methods of AVO inversion, the prestack seismic data are directly inverted into the elastic properties of the rock and an excellent lithology and fluid indicator (V_{P}/V_{S}) are provided. Then, this indicator is tested on one of the oilfields of the Persian Gulf. Moreover, by means of this method, one can locate the fluids contact and the lithological interlayers; also, by the inversion results, which are the cubes of the seismic properties of the rock, one can generate sections of the elastic properties of the rock such as Poisson’s ratio and Young modulus which are useful for geomechanical analysis. Therefore, this kind of method is a quick way for the prior analysis of the studied area. Keywords: AVO Inversion, Simultaneous Inversion, Elastic Properties, Reservoir Characterization 
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 petroelastic models (PEM) or/and fluid and lithology detection. To this end, amplitudevariationwithoffset (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.
Poststack inversion methods are the first approaches for inverting the seismic data into elastics properties. The early methods of poststack inversion were recursive methods that suffered from recovering the lowfrequency component of the reflectivity. After that, many inversion methods such as colored inversion method with spectral bluing (Neep, 2007) and modelbased 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 anglelimited 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 (poststack 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 Pwave 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 Pwave reflectivity (R_{P}), zero offset Swave 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 RauchDavies (2006) introduced more AVO attributes which are known as LambdaMuRho (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 Swave), V_{P}/V_{S}, LambdaRho, and MuRho. 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 (prestack 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 Preflectivity, Sreflectivity, and density regarding the following three assumptions:
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: Preflectivity, Sreflectivity, 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 Pimpedance, Simpedance, 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 prestack 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 poststack seismic data and prestack data. The poststack techniques transform the data to density and compressional wave velocity, which is known as Pimpedance 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 Pimpedance is affected by many parameters such as lithology, porosity, and pore fluid; moreover, it is difficult to separate their effects from each other. Although poststack 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 modelbased inversion and simultaneous inversion result are shown in Figure 1, which shows the limitation of modelbased 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 Pimpedance, Simpedance, and density are estimated, is performed (Russell and Hampson, 2006). This type of inversion falls in the second group, which is known as prestack inversion method.
a) 

b) 
Figure 1
a) Modelbased inversion result and b) simultaneous inversion. The color bars show Pimpedance 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 prestack seismic inversion are to estimates Pimpedance, Simpedance, and density, followed by predicting reservoir properties such as the property of fluid and rock. In other words, prestack inversion can be considered as a quantitative extension of AVO (Russell and Hampson, 2006). Hence, the approximated relations of Zoeppritz are needed. The reformulated Aki and Richards’ equation is given below (Fatti et al., 1994)
(1) 
where, R_{P} is zero offset Pwave reflectivity, and R_{S} represents zero offset Swave reflectivity; R_{D} is also density reflectivity.
Furthermore, the linearized zero offset Preflectivity, Sreflectivity, and density reflectivity are given by (Aki and Richards, 2002):
(2) 

(3) 

(4) 
where, , , and are difference of Pwave, Swave velocity, and density of the twolayer structure respectively, and V_{p}, V_{s}, and ρ are the averages of Pwave, Swave velocity, and density of the twolayer structure respectively. After extracting the reflectivity in Equation 1, they can be converted by poststack inversion methods such as modelbased 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):
Therefore, their algorithm uses angle gather data and allows inverting directly for Pimpedance, Simpedance, and density. This method is referred to as simultaneous inversion. They first developed it for poststack data and then extended it for prestack seismic data. The small reflectivity approximation of the Pwave reflectivity is given by (Hampson et al., 2005)
(5) 
where, i represents the i^{th} interface between layers i and i+1, and Z_{P} is Pimpedance. Thus, Equation 5 can be written in a matrix form as follows (Hampson et al., 2005):

(6) 
where, Then, if the seismic trace can be considered as a convolution of the seismic wavelet, the Earth’s reflectivity is expressed by:
(7) 
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 Pimpedance as defined by (Hampson et al., 2005):

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

(9) 
where, and Moreover, there is a linear relation between L_{S}, L_{P,} and L_{D} as follows (Hampson et al., 2005)

(10) 

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

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

(13) 
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 lowfrequency component that will not recover. Another method for solving this problem is to build an initial Pimpedance 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 OligoMiocene 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 mediumgrained 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 NWSE trending anticline, dipping by 11.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 NWSE 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 Pwave velocity, density, and gamma ray. Only one of them has Swave velocity measurement, so the estimate of Swave velocity is required for the rest of them. Therefore, the Swave velocity was calculated by means of Xu and White’s (1995) relation. In fact, the estimates of Swave velocity may be obtained from known mineral matrix properties and measured porosity and clay content or from measured Pwave 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 Pwave velocity logs. The processing steps of seismic prestack 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, Pwave, Swave, density, total porosity, and shale volume.
a) 
b) 
Figure 4
Comparison between a) conventional AVO gradient analysis and b) true amplitude DMO.
Table 1
The processing scheme of the seismic data.
SEGD read and data editing 
Navigation merge 
Deterministic zero phasing 
Lowcut filter 4 Hz 
Swell noise and interference attenuation 
Forward τ−p transformation 
Predicted deconvolution 48/180ms gap/operator 
Inverse τ−p transformation 
Old velocity input 
NMO 
kfilter spatial resampling 
CMPs to 75 m trace spacing 
f−x trace interpolation 
Static correction 
DMO 
Prestack time migration with minimum velocity 
Inverse NMO with old velocities 
New velocity picking 
NMO with new velocities 
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 40fold prestack 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 qualitycontrolling 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 Simpedance and density. Alternatively, the velocities of the rock for a more accurate analysis were inverted (Figure 8a). Therefore, it is suitable to plot the crossplot of Pimpedance with respect to V_{P}/V_{S}to determine lithology and fluid indicators. By means of these inverted elastic properties, one can also convert them to LambdaMuRho (LMR) parameters for better fluid detection. In the LMR analysis, the LambdaRho is a good attribute for this purpose. This fluid detection for sandstone member is illustrated in Figure 8b.
a) 

b) 

c) 
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) crossplot of ZS versus ZP and b) crossplot 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.
a) 

b) 
Figure 7
Inversion results; a) inverted Simpedance 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.
a) 

b) 
Figure 8
Fluids contact detection; a) V_{p}/V_{s} ratio and b) LambdaRho; 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 PP reflection prestack 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 oilwater 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 Pimpedance. From the viewpoint of detection ability, the Pimpedance, on its own, is not a gas sand or in general hydrocarbon indicator. However, by combining the ratio of the inverted P and Simpedances, which gives the V_{P}/V_{S} ratio, the interpretation becomes clearer and more efficient. Along with a drop in Pimpedance, 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, Swave 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 LambdaRho, this fluid contact is obvious (Figure 8b). As expected, the differences between V_{P}/V_{S}and LambdaRho is a key for fluid contact detection. A dolomitecemented 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 crossplot of the Pimpedance 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 hydrocarbonbearing sand has a lower Pimpedance 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.
a) 

b) 

Figure 9
a) crossplot of velocities ratio versus Pimpedance for one of the wells in the studied area in which the red, yellow, and green colors show hydrocarbonbearing sand, sandstone, and dolomite respectively; b) selected zones from a) are shown which are consistent with geological report.
a) 

b) 
Figure 10
a) crossplot of V_{P}/V_{S} versus Pimpedance, 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).
a) 

b) 
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 crossline 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 wavepropagation effects known as data uncertainties. Additionally, several uncertainties are associated with the inversion procedure such as the nonuniqueness of the inversion, which means there are several solutions which can provide the same seismic response; the need for a prior lowfrequency 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 Pimpedance, 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 prestack 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 prestack 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
AVO : Amplitude variation with offset PEM : Petro Elastic Model EEI : Extended elastic impedance : PP reflection coefficient A : AVO intercept B : AVO gradient NI : Normal incidence PR : Poisson’s reflectivity L_{P} : Natural logarithm of Pimpedance L_{S} : Natural logarithm of Simpedance L_{D} : Natural logarithm of density Greeks : Difference : Density : Incident angle (degree) : Shear modulus : Incompressibility modulus : Velocities ratio Subscripts i : Number of layer N : Number of sample 