Seismic-fluid detection-a review

Since the advent of seismic imaging techniques, the dream of geophysicists has been to identify the fluid effect and be able to accurately map hydrocarbon from the brine within a target reservoir. The usage of bright spots (strong reflection amplitudes) as an indicator of hydrocarbon was the earliest recognition of the direct role played by the pore fluids in seismic signatures. Further development of new techniques had a strong correlation with the increase in computing power and advances in seismic acquisition and processing techniques. In this review, we touch upon the relevant theory developed more than 100 years ago, and then review the methods developed over five decades leading to the quantitative interpretation of seismic data for fluid detection. We also carried out a case study to compare selected fluid identification methods applied to a complex reservoir within an oil and gas field in the Barents Sea. The impedance-based methods “CPEI-Curved Pseudo-elastic Impedance” and “LMR-Lambda-Mu-Rho” inversion provided better results compared to other techniques, highlighting the critical influence anomalous lithologies have on such screening attributes.


Introduction
There is a popular saying in the oil industry "we don't drill for sands -we drill for oil". The conventional seismic exploration, however, had been an indirect way that required mapping of potential reservoir traps, which might, or might not contain hydrocarbons. The attempt to reduce this uncertainty of presence of hydrocarbon was facilitated by the introduction of digital technology and DHI (direct hydrocarbon indicators) attributes, resulting in stronger seismic amplitude as a function of the presence of hydrocarbon (Hilterman 2001;Chopra and Marfurt 2005;Brown 2010). The DHI comprised of bright spot, dim spot, amplitude reversal, and flat spot observations in seismic sections (Fig. 1).
The relevant equations and theories (e.g., Knott 1899;Zoeppritz 1919;Gassmann 1951;Biot 1956;Biot and Willis 1957) conceived well before commercial seismic techniques were developed, are central in the fluid identification procedures. The relationship between theory and technology is essential today in the field of seismic reservoir characterization with the advancements of computing power. Today various methods developed for effective elastic moduli calculations such as Voigt (1928), Reuss (1929), Mindlin (1949), Wood (1955), Hill (1963), Hashin and Shtrikman (1963) (see Fig. 2 for the timeline) are used individually, or in combination -in the form of rock physics models. Some relevant models include the friable sand model, contact

Direct Hydrocarbon Indication (DHI)
During the late 1960s, some companies operating in the Gulf of Mexico were unusually active in license bidding, which suggested proprietary or 'extra' technology (Anstey 2005;Chopra and Marfurt 2005). Later publications showed this advantage came from bright spot technology: if one used less amplitude gain correction (AGC), and displayed data at a reduced gain, one could detect zones of substantial amplitude that could be linked to the presence of gas (Tegland 1973;Brown 2011). The problem with the display at small gain lost was that most of the structural detail from the display was missed. Therefore, companies processed two sections for each play -the bright spot section and the conventional section. This had the disadvantage that, in the interpretation of bright spots, the criterion of geological credibility was made less direct, and therefore weakened (Anstey 2005).
The technique of digital recording was a milestone in processing and enabled the technologies of preserving relative seismic amplitudes. Early successes led to the rapid development of technology in the early 1970s, including efforts to quantify seismic amplitude changes and to calculate pay-sand thicknesses (Forrest 2000;Chopra and Marfurt 2005). DHI also comprised phase reversal, flat spots, frequency loss, time sags, time shadows, polarity reversals, and dim spots features identified by the interpreter that became the motivation for later seismic-attribute developments (Sheriff and Geldart 1995;Chopra and Marfurt 2005;Brown 2011). (Bortfeld 1961) Bortfeld (1961) developed approximations to the Zoeppritz (1919) equations for calculating P-and S-wave reflection and transmission coefficients by assuming small contrasts between layer properties (Eq. 1). Though Zoeppritz's equations were the exact solutions, they were complicated to fully calculate and not physically intuitive, whereas the Bortfeld (1961) equation manifested the influence of fluid and the rigid part on the reflectivity:

Approximation to Zoeppritz Equations
where, R P (θ 1 ) is reflection coefficient at an angle of incidence θ 1 , θ 2 is transmission angle, V P1 is P -wave velocity of layer 1, V P2 is P -wave velocity of layer 2, V S1 is S -wave velocity of layer 1, V S2 is S -wave velocity of layer 2, ρ 1 is density of layer 1, and ρ 2 is density of layer 2 (Fig. 3a). Bortfeld emphasized that all converted waves depended primarily on the Poisson's ratio, as Koefoed (1955) had indicated earlier. However, in areas where the interval velocities are not well known, the usage of these equations may not be accurate.

Aki and Richards (1980)Approximations
In the 1980s, Aki and Richards presented the amplitudes in a seismic gather as a linearized form of the Zoeppritz (1919) amplitude  Løseth et al. 2009). Note also the increased amplitude and continuity in the caprock just above the hydrocarbon-filled reservoir; b) Compaction trends of sand and shale in terms of acoustic impedance to show the nature of DHI as a function of depth/age (modified after Brown 2010). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (1) M. Fawad, et al. Earth-Science Reviews 210 (2020) 103347 versus offset (AVO) equations, which can be further approximated (Shuey, 1985):  This formulation added two extra terms to the zero-offset case, a gradient term G and a curvature term C, often referred to as A, B, and C, where the term A is called the intercept (Fig. 3b). M. Fawad, et al. Earth-Science Reviews 210 (2020) 103347 3.3. Ostrander (1984)AVO method Ostrander (1984) suggested a method of analysis using seismic reflection amplitude versus shot-to-group offset distinguishing between gas-related amplitude anomalies and other types of amplitude anomalies. He also recognized and emphasized that Poisson's ratio (σ) had a substantial influence on the changes in reflection coefficients as a function of the angle of incidence. Poisson's ratio is defined as: He showed that the gas sands have a lower Poisson's ratio compared to that of brine sands and that this difference can be used to identify the  Aki-Richards (1980) approximation to the Zoeppritz equations (modified after Simm and Bacon, 2014); (c) plot of P-wave reflection coefficients vs. Angle of incidence for a reduction of Poisson's ratios across an interface (modified after Ostrander 1984); (d) approximations in an instance with the velocity difference ΔV < 0, the elastic properties correspond to actual gas-bearing sand in the Gulf of Mexico (modified after Shuey 1985); (e) Three classes of gas sands based on their AVO characteristics (modified after Rutherford and Williams 1989). gas sands employing the difference of amplitude variation with offset (AVO) signatures from the two different pore phase cases (Fig. 3c). He, however, advised considering other factors that affected observed reflection amplitude versus offset and underscored the importance of amplitude balancing during seismic processing. (Shuey 1985) Shuey (1985) simplified the compressional wave reflection with offset coefficient R P (θ) given by the Zoeppritz equations (Fig. 3d), leading to the possibility of extracting density data from wide-angle AVO. Shuey's (1985) three-term approximation to the Zoeppritz equation for P-wave (V P ) reflectivity was:

AVO approximation for wide angle
where R 0 is the amplitude at normal incidence (intercept), . A 0 specifies the normal, gradual decrease of amplitude with offset. The density (ρ) contrast at an interface can be calculated by subtracting the 3rd coefficient from the intercept (Avseth et al., 2005): Theoretically, by quantifying the density, it was possible to discriminate residual gas from commercial gas. However various authors (e.g. Swan et al. 1993;Chen et al. 2001;Hilterman 2001;Avseth et al. 2005) underscored the difficulty of calculating density from three-term AVO because of the poor signal-to-noise ratio of the 3rd coefficient C, acquisition effects, processing complications with increasing offset, and anisotropy. (Rutherford and Williams 1989) Based on AVO characteristics, Rutherford and Williams (1989) grouped the gas sands into three classes defined in terms of the zeroangle reflection coefficient (R 0 ) at the top of the gas sand. Class 1 sands were defined as having higher impedance than the overlying shale with relatively large positive values for R 0 . Class 2/2P sands were having nearly the same impedance as the overlying shale and were characterized by values of R 0 near zero. Class 3 sands were considered to be having a lower impedance than the overlying shales with large negative values for R 0 . Each of these sand classes was explained to have a distinct AVO signature (Fig. 3e). Class 4 sand was identified later (Castagna and Swan 1997), and characterized by similar R 0 as that of Class 3, but the negative amplitude decreased with increasing offset. (Smith and Gidlow 1987) Smith and Gidlow (1987) came forward with the concept of "weighted stack" in which a "difference stack" relative to a wet background trend is generated.

Fluid factor
Where v is the background V S /V P ratio which can be predicted by application of the mudrock line to interval velocities obtained from conventional velocity analysis. It became possible to display the presence of gas on a stack without any post-stack processing using this procedure (Fig. 4a).

Fluid
Factor in term of V P and V S (Fatti et al. 1994), and Reflection difference (Castagna and Smith 1994) The reflection difference (Castagna and Smith 1994) and the fluid factor trace (Fatti et al. 1994) were both constructed so that all reflectors associated with the rocks deviating from the mudrock line (Castagna et al. 1985) would show bright amplitude. Fatti et al. (1994) defined the fluid factor in terms of P-and S-wave reflectivity: Furthermore, the "F-impedance" traces ( Fig. 4b) can be obtained by: where X and M are the intercept and the slope of the mudrock line, respectively.
According to Castagna and Smith (1994), the AVO product indicator (A*B) was often inappropriate for hydrocarbon detection. They proposed that in clastic stratigraphic sections, the normal P-wave and Swave reflection coefficient difference R P -R S is applicable when corrected for overburden effects. They put forward the following equation: where A = R 0 , the normal incidence P-wave reflection coefficient, and the AVO gradient (slope) can be expressed as: The R P -R S is always negative for shale over gas-sand interfaces and significantly more negative compared to the shale over brine-sand interfaces (Fig. 4c). The R P -R S works for any gas-sand impedance providing there is a strong gas effect on the P-wave velocity. However, in both procedures (Castagna and Smith 1994;Fatti et al. 1994) several lithologies, which do not follow the Mudrock line, will also brighten up on these attributes. (Verm and Hilterman 1995) Verm and Hilterman (1995) carried out a further approximation of Shuey's (1985) equation resulting in:

Poisson reflectivity
where NI = R 0 , which is the normal incidence reflectivity, and PR is the Poisson reflectivity = Δσ/(1-σ) 2 . By colour-coding the NI PR matrix, discrimination of fluid bearing lithologies became possible as the separation of reflection clusters increased (Fig. 5a). Furthermore, the ability to carry out a 45° rotation of the NI and PR improved the delineation of low contrast reservoir interfaces.

Intercept versus Gradient cross-plot (Castagna and Swan 1997)
Castagna and Swan (1997) suggested that the hydrocarbon-bearing sands should be classified according to their location in the AVO intercept (A) versus gradient (B) cross-plot. A background trend within the zone of interest must be defined with the help of well control if the seismic data are correctly amplitude calibrated or with the seismic data excluding hidden hydrocarbon-bearing zones. Upon plotting, the top of gas-bearing sands plot under the background trend, whereas the bottom of gas sand reflections plot above the trend. Their classification was identical to that of Rutherford and Williams (1989); however, they subdivided the Class 3 sands (low impedance) into two classes (3 and 4). The Class 3 and 4 gas sands may have identical normal incidence reflection coefficients, but the magnitude of Class 4 sand reflection coefficient decreases with offset or increasing angle of incidence while Class 3 sand reflection coefficient increases. The dimensionality reduction benefit of this qualitative method is that the information on cross plots can be reduced to one-dimensional parameters, which are easier to interpret (Fig. 5b). The residual gas can give false anomalies, which may become difficult to differentiate from that of commercial gas. A problem encountered with this approach is that the background, or the lithology trend, is usually much steeper than what the theory based on rock properties predicts. According to Chopra and Castagna (2014), the trends in Intercept-Gradient cross-plots can be caused by lithology in some cases and other cases, by noise. According to de Bruin (2019) the very strong apparent trend is created because while displaying amplitudes in the Intercept-Gradient plane, we ourselves (rather than lithology or noise) unintentionally decide how steep the trend will be. This decision depends on the function we chose to which we match the prestack amplitudes. To make sure if a meaningful trend exists in a prestack data he recommended to plot a near and far amplitude cross-plot to check for two data clusters. (Goodway et al. 1997) Another form of Fatti et al. (1994) equation relating the pre-stack seismic data with the acoustic impedance reflectivity (ΔAI/2AI), shear impedance reflectivity (ΔSI/2SI) and density reflectivity (Δρ/2ρ) term is:

LMR (Lambda-Mu-Rho) inversion
where AI = ρV P and SI = ρV S . For the angle of incidence (θ) less than 35° and V P /V S ratio between 1.5 and 2.0 the third term in the above equation can be neglected (Gidlow et al., 1992).
For LMR inversion (Goodway et al. 1997) volumes of zero-offset Pwave and S-wave reflectivity are estimated through least square fitting of R P to the above equation. Volumes of P-Impedance (AI) and S-Impedance (SI) are then obtained through blocky model-based inversion, constrained by horizons, sonic logs, and an initial guess model (Young and Tatham 2007) (Fig. 5c).
The λρ and μρ volumes are obtained using the following relationships: These attributes are cross plotted (Fig. 5c), and zone analysis is performed to identify a zone of interest in cross-section or map view. Since the impedances are estimated from seismic inversions, there is a possibility of error introduction, furthermore squaring the impedances and taking their linear combinations introduces further errors and bias in the Lambda-Rho and Mu-Rho estimates (Avseth et al. 2005). To avoid this, (Gray et al (1999)) and (Gray (2002)) suggested using an approximation for R P expressed directly in terms of λ and μ.

Acoustic impedance inversion (Lindseth 1972)
Lindseth (1972) developed a new method of processing, displaying, and interpreting seismic data, which increased the amount of geologic information (i.e., lateral changes in lithology and porosity) available to the interpreter. Seismic field data were deconvolved by wavelet processing, inverted and combined with low-frequency velocity components derived from reflection analyses to produce synthetic sonic logs. The procedure initially generated an inverted impedance section that Lindseth (1972) converted to a sonic log section using a general empirical relation. The basic procedure of inversion of seismic trace data to generate impedance logs was first described by Delas et al. (1970), published by Lindseth (1972), and Lavergne (1975) etc. This was mainly the inverse of producing seismograms from borehole sonic logs  Castagna and Swan 1997); (c) LambdaRho versus MuRho crossplot for discriminating gas-bearing sands from brine sands, shales and carbonates (modified after Goodway et al. 1997). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) (Lindseth 1979) (Fig. 6a). The method received wide acceptance for being prompt and accurate, and was used in stratigraphic interpretation framework that picked up at that time (Chopra and Marfurt, 2005).
In a stacked seismic section the traces can be modelled as the convolution of the earth's reflectivity and a bandlimited seismic wavelet (Russell and Hampson 2006), which can be written  where S t is the seismic trace, W t is the seismic wavelet, r t is the reflectivity and '*' denotes convolution. The zero offset reflectivity, in turn, is related to the acoustic impedance of the earth by where R 0i is the zero-offset P -wave reflection coefficient at the ith interface of a stack of N layers and AI i = ρ i V Pi is the ith P-impedance of the ith layer, where ρ is density and V P is P-wave velocity. Assuming that the eq. (16) represents the recorded seismic signal (Lindseth 1979), this can be inverted using the recursive equation to recover the P-impedance/acoustic-impedance (AI).
So, by applying eq. (17) to a seismic trace we can effectively transform, or invert, the seismic reflection data to P-impedance. Lindseth (1979), however, recognized some problems with this procedure (Russell and Hampson 2006) such that the recorded seismic trace is not the reflectivity shown in eq. (16) but rather the convolutional model given in eq. (15). Furthermore, the effect of the bandlimited wavelet is to remove the low-frequency component of the reflectivity that cannot be recovered by the recursive inversion procedure of eq. (17).
The original recursive or trace-integration seismic inversion technique for acoustic impedance later evolved during the late 1980s and 1990s, which progressed to model-based inversion, sparse-spike inversion, stratigraphic inversion, and geostatistical inversion providing reasonably accurate results (Chopra and Kuhn 2001).

Elastic Impedance
Connolly (1999) introduced elastic impedance (EI), which computes conventional acoustic impedance for a non-normal angle of incidence. EI allowed relevant well log data to be tied directly to the far offset seismic observations, which could then be calibrated and inverted without reference to the near offsets. It was defined as: where K is a constant, set to the average of (V S /V P ) 2 over the interval of interest. Fig. 6b shows the difference between AI and EI (30°) from well log data in a discovery well and core sample measurements in a well in the Foinanven Field (Connolly 1999). It is clear that the far offsets are more sensitive to changing saturation than the near offsets. Using this procedure, the first order AVO effects can be incorporated in seismic enabling extraction of quantitative AVO information from large 3D volumes. The limitation is that as |Sin 2 θ| approaches 1, the EI log becomes increasingly inaccurate. This problems was overcome by introducing constants V Po , V So and ρ o within the EI function (Whitcombe 2002) : The modifications allowed a direct comparison between elastic impedance values across a range of offset angles.  Using an extension to the elastic impedance concept, Whitcombe et al. (2002) obtained the optimum projection for a noise-free, isotropic environment. The elastic impedance is itself an extension of acoustic impedance (AI) to nonzero angles of incidence. To obtain extended elastic inversion, they modified the definition of elastic impedance (EI) beyond the range of physically meaningful angles by substituting tanχ for sin 2 θ in the two-term reflectivity equation:

EEI (Extended Elastic Inversion)
The primary variable now becomes χ rather than θ. K is a constant, normally set to the average of (V S /V P ) 2 over the interval of interest. The normalizing constants V P0 , V S0 , and ρ 0 can be taken to be either the average values of velocities and densities over the zone of interest or the values at the top of the target zone. With this normalization, the elastic impedance has the same dimensionality as the acoustic impedance. The EEI is defined over angle χ (chi) ranging from −90° to +90°. It should not be interpreted as the actual reflection angle but rather as the independent input variable in the definition of EEI. That enabled the identification of different areas of EEI space optimum for fluid and lithology imaging (Fig. 6c). Once the appropriate χ value is obtained, the equivalent seismic section can be generated while routine AVO processing from combinations of intercept and gradient stacks.

Simultaneous Inversion (Ma 2001)
Ma (2001) combined the AVO extraction and impedance inversion into a single step, and transformed P-wave offset seismic gathers into Pand S-impedances using simulated annealing. Instead of reflection coefficients the use of impedances as model parameters, allowed reliable and flexible constraints to be included on the inversion algorithm.
The Fatti et al. (1994) equation can be written as: Ma (2001Ma ( , 2002 replaced the average Vs/V P ratio by SI/AI ratio, so that the reflection coefficient R(θ) became only related to three parameters, AI, SI and θ. The angle of incidence (θ) could be calculated using a ray-tracing method (Smith and Gidlow 1987). This assumption is good when the density change between the two adjacent layers is small. The SI/AI ratio was not determined from background P-and S-Impedance model, but was estimated from the impedance models at each iteration. So, based on the following equation a simultaneous inversion procedure was proposed to derive acoustic and shear impedances from prestack offset seismic gathers: The assumptions were that the earth contains approximately horizontal layers at each common depth point (CDP), and both acoustic and shear impedances describe each layer. The simultaneous inversion was achieved by using the Monte-Carlo approach in the form of simulated annealing. An example of results of simultaneous inversion is shown in Fig. 7a. Hampson et al. (2005) modified the Fatti et al. (1994) formulation using a small reflectivity approximation coming up with the following equation: where T(θ) represents a seismic trace at a given angle θ, while W and D represent the angle-dependent wavelet and the derivative matrix respectively. The terms Δln(SI) and Δln(ρ) represent a deviation from the background linear trends. The term ã 1 is modified from the original a 1 = 1 + tan 2 θ using regression coefficients from the background trend,

Bayesian linearized AVO inversion (Buland and Omre 2003)
Due to the commonly highly variable nature of the various reservoir properties contrasting to, but often being represented by consistent average behaviour, it was imperative to combine deterministic physical models with statistical techniques. This approach led to new methods for interpretation and estimation of reservoir rock properties from seismic data. The procedure identified the most likely interpretation, the uncertainty of the interpretation, and then guided the quantitative decision analysis (Mukerji et al. 2001;Avseth et al. 2005). Buland and Omre (2003) defined as inversion in a Bayesian setting where the complete solution was a posterior model. The prior model is expressed mathematically by a probability density function (PDF). The posterior distributions for P-wave velocity, S-wave velocity, Density, Acoustic impedance, Shear impedance and V P /V S ratio can be obtained (Fig. 7b). The probability for lithology and fluid classes can be obtained in a similar way. The inversion algorithm is based on the convolutional model and a linearized weak contrast approximation of the Zoeppritz equation (Aki and Richards 1980). Acoustic impedance is the best determined parameter; however, the inversion was not good at getting the density information. Generally, the resolution of the various parameters was a function of the prior model and the noise covariance. (Johansen et al. 2013) Seismic reservoir characterization is the conversion and combination of seismically derived elastic properties such as P-and S-wave velocities, acoustic impedances, elastic impedances, or other seismic attributes into properties related to lithology and reservoir conditions. A variety of rock physics models exist to find these relationships. The rock physics models, however, are specific for the type of lithology, porosity range and shape, texture, saturation, and the pore fluid conditions. There are usually too many rock physics parameters compared to the seismic parameters; this is known to be an underdetermined problem with non-unique solutions (Johansen et al. 2013). The procedure of inverse rock physics modelling which aimed at direct quantitative prediction of lithology and reservoir quality from seismic parameters, with parallel handling of non-uniqueness and data error propagation. It is based on a numerical reformulation of rock physics models so that the seismic parameters are input and the reservoir quality data are output. The modelling procedure can be used to evaluate the validity of various rock physics models for a given data set. Furthermore, for a given rock physics model, it delivers the most robust data parameter combinations to use for either porosity, lithology, or pore fluid prediction (Johansen et al. 2013).

Inverse Rock Physics Modelling
Inverse rock physics modelling relates how well the input data match the rock-physics model for a particular prediction of rock properties. The search is done in forward-modelled so-called constraint cubes (Johansen et al. 2013;Bredesen et al. 2015) relating rock parameters to seismic parameters (Fig. 8a). The strength of this procedure is the estimation of the probability distribution of the predicted parameters. However, predictions in frontier basins may be uncertain outside the range of the well control. (Avseth et al. 2014) CPEI is pseudo-elastic impedance that mimics the EEI using the AI and V P /Vs directly. As the target zones are located at different burial depths over a larger area, the estimate of values from a nonlinear curve representing the rock-physics model itself is suggested. This nonlinear attribute is defined as:
where, AI T is the tangent of the rock-physics model at a point where AI and V P /V S represent mean values for water saturated sandstones. V P0 / V S0 represents the intercept of the function with the V P /V S axis, whereas a, b, c and d are the constants controlling the slope of the line. Also, γ is a scaling parameter that makes the PEI values having the same range as EEI at chi angle χ = 24.
The CPEI attribute derived from AI and V P /Vs inversion data is shown in Fig. 8b, the random line intersects two wells. The easternmost well (Well 1), a gas-discovery, is located on a structural high, and the gas reservoir in Middle Jurassic is clearly seen in this attribute. Well 2 encountered a thick Upper Jurassic hot shale and only a thin watersaturated Middle Jurassic sandstone. The CPEI attribute shows no anomaly at this well (Avseth et al. 2014).
This algorithm was suggested to screen seismic data for the identification of hydrocarbon leads. This equation can be used for a Fig. 8. (a) Forward-modelled rock-physics constraint cubes for AI (left), V P /V S (middle) and corresponding solutions for AI = 6 km/s × g/cm3 and V P /V S = 1.6 (right) with the colour scale indicating model probability (modified after Bredesen et al. 2015); (b) random line of CPEI optimized for fluid effects intersecting Well 1 and Well 2. A strong fluid anomaly is evident corresponding to the proven gas reservoir in Well 1 below the (blue) Base Cretaceous Unconformity (BCU) horizon. The inserted cross-plot is derived from a 200-ms interval beneath the BCU horizon (Avseth et al. 2014); (c) Comparison of the upscaled actual and predicted S w log in the √Is -1/σ plane for Well-A. (left), comparison between the upscaled actual and predicted S w logs for the same well (middle), a resultant S w volume cross-section along Well-A and Well-B (blind well) together with its S w log (right). There is a good match between the seismic and well-log derived S w (modified after Alvarez et al. 2015); (d) density ratio map (upper panel) obtained from Zoeppritz's inversion and the corresponding AI r -ρ r scatterplot (blue points in the lower panel). The critical angle contours are calculated from eq. 27. The anomaly (approximately) resides in the ρ r/ AI r ≥ 1 domain; furthermore, most of the normalised kernel density estimate (NKDE) (Lehocki et al. 2015) points just plot around the intersection of the AI r = 1 and ρ r = 1 lines (modified after Lehocki et al. 2019). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) qualitative indication of hydrocarbon, however the parameter adjustment to fit corresponding wet sandstone trend and reservoir matrix mineralogy is cumbersome.

Multi-attribute rotation scheme (MARS) (Alvarez et al. 2015)
The multi-attribute rotation scheme (Alvarez et al. 2015) used a numerical solution to estimate a transform to estimate petrophysical properties from elastic attributes. This was attained by estimating a new attribute in the direction of maximum change of a target property in an n-dimensional Euclidian space formed by n attributes, and subsequent scaling of this attribute to the target unit properties. This approach was carried out using well log derived elastic attributes and petrophysical properties and applied over seismically-derived attributes to generate volumes of those properties. For the specific case of two dimensions the equation used was: where, A1 and A2 are elastic attributes; SF A1 , SF A2 are scale factors, which are applied to equalize the order of magnitude of the attributes; and θ i is the angle where the maximum correlation is reached. As different attribute combinations produce different results, it is possible screening n-dimensional spaces formed by n-number of attributes and angle evaluations for each one of these spaces, with the objective of finding an attribute τ that represents the global maximum correlation with the target petrophysical property.
For a case study, MARS was used to estimate a water saturation (S w ) volume in a mud-rich turbidite gas reservoir located in onshore Colombia (Alvarez et al. 2015). Among the two wells (Well A and Well B), Well A was used for MARS analysis. The other well was used as blind test of the technique. Finally, on applying the resultant transform to seismically-derived volumes of √Is and 1/σ to obtain a volume of S w . The well log analysis and the resultant S w volume cross-section along the two wells with their S w logs is shown in Fig. 8c. Lehocki et al. (2019) suggested an inversion of Zoeppritz equation (Zoeppritz 1919) to obtain the ratio of density of two layers at the layers' interface. The input parameters to the PP-Zoeppritz equation for plane waves can be written as: , index 1 refers to layer 1, index 2 refers to layer 2, and θ is the angle of incidence of the plane wave.

Density ratio inversion from Zoeppritz equation (Lehocki et al. 2019)
The critical angle θ c that occurs at which the refracted wave moves along the interface between the two layers is represented by (Lehocki et al. 2015(Lehocki et al. , 2019: The PP-Zoeppritz equation (eq. 26) can be inverted to extract the density ratio: Since the Zoeppritz equations are highly nonlinear, and it is challenging to excavate the density ratio, the inversion leads to a 12th degree polynomial equation. With today's computational power such level of calculations are manageable.
The discrimination between a low saturation gas reservoir (fizz gas) and commercial gas is often difficult using the existing fluid identification techniques. The distinction seems possible (Fig. 8d) employing the density ratio technique even in (initially) cemented rocks (Lehocki et al. 2019) as the diagenetic cement dampens the fluid effect on elastic properties.

CASE STUDY -A Barents Sea example
We attempted to test some amplitude/reflectivity derived and impedance-based fluid identification techniques on a complex reservoir in the Goliat oil and gas field to compare their results. The Goliat field is a faulted structural closure in the crestal part of a major northeastsouthwest trending roll-over anticline located in the southeastern part of the Hammerfest Basin within the Norwegian Barents sea (Yenwongfai et al. 2017a(Yenwongfai et al. & b, 2018Mulrooney et al. 2018; NPD, 2019) (Fig. 9a).
A high-quality seismic dataset courtesy of Vår Energi AS and license partners covering the Goliat field was utilised to test and compare selected methods for seismic fluid characterization. These methods include both AVO attributes (amplitude) and elastic attributes derived from deterministic inversion for P/S-impedance. We employed seven wells, multi-azimuth (MAZ) 3D seismic data and prestack depth-migrated velocities. Angle stacks corresponding to 17° (near), 32° (mid) and 45° (far) formed the basis of AVO analysis and served as input for prestack elastic inversion. The data were processed with consideration of AVO workflows (Buia et al. 2010), and are thus suitable for testing seismic fluid characterization. We employed the Hampson-Russell software package to process data and perform calculations using the workflows associated with the fluid detection methods described in the previous section.
Upper Triassic-Jurassic reservoir intervals of the Realgrunnen Subgroup (Fruholmen and/or Tubåen formations) contain oil and gas in the Goliat field, and are capped by Upper Jurassic organic-rich shales of the Fuglen and Hekkingen formations (Fig. 9b). In this well (No. 3), the Fuglen Formation is only 14 m thick, with V P , V S and density higher than that of the underlying reservoir. The Fuglen Formation's thickness close to the seismic resolution and high magnitude of elastic properties have potential to influence the seismic signature at the reservoir-seal interface. A present reservoir burial depth of around 1 km with an approximate net Cenozoic uplift of 1-1.2 km means that the sandstones encountered are expected to be consolidated and at least modestly quartz cemented unless clay-coated. This can negatively affect the fluid sensitivity, but the calculated porosity in clean sandstone intervals still approaches 30-35% in some cases.
Comparing the oil-and gas-bearing zones within the reservoir (Fig. 9a), it is evident that in the southwestern part of the field this interface constitutes shale over gas sand (well 3), whereas in eastern and northern parts of the field, the reservoir is oil-bearing (wells 1, 2 and 6). A recently drilled well (well 7) also identified a gas cap in a disconnected fault segment to the northwest, but our information from this wellbore is limited. Well 5 is dry, whereas well 4 on the southeastern side is water-bearing with oil shows within the Realgrunnen Subgroup.

Results
For representation, the amplitude/reflectivity derived methods of fluid indication were extracted along the top reservoir reflector. The inversion based attributes were computed as mean value over a 20 m window from the top reservoir reflection and are represented in map view. All methods are exercised following standard workflows and procedures. Fig. 10a-c displays a comparison of different expressions of angledependent amplitude in map view. A first-order observation is that the reservoir is represented by dimming of amplitudes, for instance visualised by the near*far reflection amplitude (Fig. 10a), near-far reflection amplitude (Fig. 10b), and Intercept*gradient (Fig. 10c). The top reservoir reflection AVO signature shown in the inserted θ versus R P (θ) crossplot (Fig. 10a) indicates that class IV AVO is observed around the gas sand (with phase change at high angles; near well 3, as well as outside the reservoir (near well 5). The top oil-bearing reservoir reflector, however, classifies as weak AVO class II (near well 1). It is obvious that near/ far stack operations and intercept-gradient product do not yield meaningful results owing to predominantly class IV AVO.

AVO Mapping
The fluid factor (Fig. 10d) attribute performs poorly in highlighting the extent of the oil-bearing part of the reservoir. The highest fluid factor is observed from the area around well 3 with gas and toward southwest, with similar warm colors in the down-faulted area near the fault relay intersection and westward (around well 5). Seeing as the attribute was developed to identify gas, this is no surprise, as the softening due to low sand content (around well 5) appears to be greater than the marginal softening associated with oil saturation in cemented sand (wells 1, 2, and 6).

Seismic Inversion
The main observable difference when instead considering a window below the top reservoir reflection (seismic inversion attributes), is that the gas-bearing region of the reservoir is more.
clearly highlighted, where a suitable window size is chosen (Fig. 11a-e). Properties derived from impedance inversion include Pimpedance (Fig. 11a) and S-impedance, which transfers to V P /V S ratio (Fig. 11b) and Poisson's ratio (Fig. 11c). We can note meaningful correlation between all three properties and the field outline, but similar values are also detected outside the main reservoir. Evidence of this is observed in the calibration point in well 5. Note also that values near the gas-filled section is characterized by similar values as regions outside the rest of the reservoir.
Extended elastic impedance inversion (Fig. 11d) based on the EEI correlation analysis of Yenwongfai et al. (2017a) displays a clearly different reservoir prediction (bright warm colour) compared to the methods shown in Figs. 10 and 11a-c. The EEI correlation coefficient varies between 0.6 and 0.9, of which the highest was found in well 3 (gas+oil), and the lowest in well 1 (oil).
The EI method (Fig. 11e-f) is unable to even locate the gas reservoir from the seismic data, even when carefully considering information from well logs. The gas anomaly is not defined by a unique V P /V S relationship compared to brine-filled lithologies, as seen in Fig. 11f with well log data from the gas well. The range of values observed in the well is not represented anywhere in the seismic, and the EI method is not advanced enough in this case since it builds on a simple difference in V P /V S from a mudrock (brine) trend.

Crossplot techniques
The λ-μ-ρ crossplot in Fig. 12a shows how the standard cutoff for fluid indication of 20 GPa*g/cm3 superimposed on well log data representing gas sandstone, oil sandstone (similar to silty/shaly reservoir sections) and organic-rich shale (well 3). We note how these properties fail to distinguish organic-rich shales from gas-bearing sandstone, as the soft source rock data falls close to gas sand data (see cross plot in Fig. 12b). We see the implied accuracy of gas sand identification from the seismic data near a few control points in cross-section view in Fig. 12b, using the LMR attribute and crossplot zonation. The match is good between gas indicators and knowledge from wells, but oil-bearing sand and silty/shaly reservoir formation intervals are too similar to separate with confidence using this generalized approach.

Recent Developments
In the CPEI attribute plot ( Fig. 12c; < 7.2 g/cm3 × km/s indicative of HC), we note that the oil-bearing part of the field blends in with the background data (brine saturation and clay-rich reservoir formations) due to the averaging window. This seems similar to the LMR attribute based on crossplot classification of fluid/lithology (Fig. 12a). Since the CPEI method relies only on the distance of points from a reference wet sandstone curve onto an AI-V P /V S plane, the gas anomalies are very prominent because of low AI and V P /V S ratio within the gas bearing sandstone. In a combined oil and gas reservoir there is however no possibility of determining a universal cutoff separating the oil and gas effects with CPEI.

Discussion
A significant issue while testing fluid identification methods on Jurassic reservoirs on the Barents Sea is the presence of organic-rich shale capping the reservoir. Low acoustic impedance and V P /V S ratio exhibited by the organic shale caprocks very often imitate the oil bearing sandstone seismic signature, or potentially a tight gas sandstone. This behaviour of the organic shales is typically very different from the grey, conventional shale/claystone assumed to constitute the dominant background lithology in older methods, from which the hydrocarbon-bearing reservoirs are meant to stand out. Internally in the caprock section, there are also large differences in elastic properties between Hekkingen and Fuglen Formation (crossplots in Figs. 11f  and12b). These features are also evident from the AVO behaviour, where the diverse caprock properties is the main deciding factor for the AVO class observed rather than reservoir fluid content. Observe the  Mu-rho crossplot in the complex reservoir-cap rock system; (c) CPEI attribute map with hydrocarbon sand cutoff of 7.2 g/cm 3 × km/s. difference between CPEI versus the AVO based indicators. In a siliciclastic sequence with absence of coal or limestone, the fluid factor is theoretically derived from the difference of two reflectivities from the top interface of a very thick reservoir, however, it is independent of the medium properties above the interface subject to ignoring the higherorder term of reflectivity (Zhou and Hilterman 2010). The typical requirement of assuming a "background" V P /V S ratio or relating the method to the mudrock line is basically the simplistic way of accounting for different compaction levels/stages. However, the more divergent are the compaction trends, or different lithologies are present, the less successful this approach will be. This explains why AI-V P / V S trends (CPEI) are better defining particularly the gas reservoir sections. The complex geology is represented by the difficulty in separating data clusters, as seen for example in the EI crossplot (Fig. 11f).
The elastic inversion attributes particularly highlight how fluid sensitivity, in the case of oil, is compromised by what is most likely a few percent of cement in the reservoir sandstone. Considering depthtrends of different lithologies and basing the fluid discrimination on corresponding rock physics templates (such as the CPEI method) results in a better segregation of what roughly corresponds to the gas-capped reservoir, but we still struggle with discrimination of oil-saturated sands. This is also reflected in low correlation values between EEI and saturation for wells only encountering oil (Yenwongfai et al. 2017a). The EEI method likely also suffers more from being unable to decouple lithology and saturation.
We are able to use crossplot classification to make fairly good predictions of reservoir fluid content, but such approaches require good a priori knowledge about the target structure, e.g., as provided by wells. It is clear that older methods that are based on separating highly contrasting gas sands from grey shales are not advanced enough to work independently in complex geological settings such as that observed on the Barents Sea without considerable information about the reservoir in advance. This observation especially holds true if we consider a "screening" type application of seismic fluid detection methods. A statistical approach could have benefits in such datasets where advanced clustering could be required to separate lithologies and fluids.

Conclusions
Over the last five decades, we have seen the fluid detection procedures ranging from qualitative bright spots to relatively quantitative cross-plotting and seismic section generation. With the progress of both computer power and advanced seismic acquisition and processing techniques, fluid identification and quantification techniques have been gradually moving more toward a quantitative approach.
The selected attributes (partial stack operations, AVO interceptgradient, Fluid Factor, CPEI, LMR, EI and EEI) tested on a Barents Sea reservoir capped by organic-rich shales exhibited varied results. The reflectivity based methods (i.e., Partial Stack Maths, AVO interceptgradient, Fluid Factor generally did not yield meaningful anomalies in that complex reservoir-seal setup possibly owing to presence of thin Fuglen organic-rich shales above the reservoir sandstones. The impedance-derived methods CPEI and LMR gave better hydrocarbon representation, however the EI and EEI did not achieve comparable results due to poor ability to decouple saturation and various lithologies, and particularly for EEI the relatively poor correlation with saturation in the case of oil. In general, the presence of organic-rich shales above the reservoir level is important to consider before attempting remote fluid detection.

Declaration of Competing Interest
None.