Poroelastic property analysis of seismic low-frequency shadows associated with gas reservoirs

Seismic low-frequency amplitude shadows have been widely used as a hydrocarbon indicator. This study investigates the effect of reservoir properties and seismic wave mode conversion on the characteristics of the low-frequency amplitude shadows in gas-bearing reservoirs. The target gas reservoirs are typically related to the lithology of tight sandstone with strong heterogeneity. Pore-fluid distribution within the reservoirs presents patchy saturation in the vertical and horizontal directions, and this patchy saturation easily induces low-frequency shadows beneath gas-bearing reservoirs. These low-frequency shadows are validated by using a poroelastic simulation method. The results of our field case-based study indicate that pore-fluid property, plus the thickness and heterogeneity of reservoirs are the key elements in the generation of low-frequency shadows. The results also indicate that the poroelastic simulation method can be used to effectively predict the spatial distribution of gas-bearing reservoirs, by directly verifying the low-frequency shadow phenomenon existing in the seismic data.


Introduction
Seismic low-frequency shadows have been widely used as a hydrocarbon indicator (Taner et al. 1979). A low-frequency shadow (LFS) is a seismic reflection with a strong amplitude at low frequencies immediately beneath a hydrocarbonbearing reservoir. Numerous application examples of LFSs as hydrocarbon indicators have been reported (Goloshubin et al. 2006;Wang 2010;Aghdam & Riahi 2015). The goal of this study was to investigate how the physical property of reservoirs and the conversion of seismic wave modes affect the characteristics of LFSs, and the investigation was based on a field case study.
A strong low-frequency amplitude anomaly zone near reservoirs cannot be interpreted by the seismic attenuation mechanism alone, which simply does not enhance the lowfrequency amplitude (Xu et al. 2011;Rao & Wang 2015). Apparently, the strong amplitudes within the low-frequency anomalies are increased by certain physical processes. Ebrom (2004) lists at least 10 causes that may lead to the generation of LFSs, such as seismic data processing-related factors, attenuation and wave interferences. Carcione et al. (2018) use numerical modelling of LFSs to quantify the effects of seismic attenuation and stretching of normal moveout (NMO) on the generation of LFSs. Pride et al. (2002) report that Biot's slow wave is likely to be the dominant source of low-frequency anomalies in a sedimentary basin, and the layer thickness significantly affects the amplitude of the slow P-wave. Chabyshova & Goloshubin (2014) note that LFSs are due to the converted slow P-wave that is generated in thinly layered permeable fluid-saturated media. Results from field data and theoretical simulations show that LFSs may be caused by multiple factors. It is indispensable to investigate the influencing factors that may alter the characteristics of the LFS, if applying the LFS analysis to field data for the purpose of identifying gas reservoirs.
The prevailing perspective is that the low-frequency anomalies are mainly generated due to the fact that highfrequency components are attenuated more severely than the corresponding low-frequency components, such as seismic waves propagating through fluid-saturated media (Korneev et al. 2004;Quintal et al. 2009;Liu et al. 2011;Dupuy & Stovas 2015;. However, mesoscopic loss effects due to partial saturation, either with gas or with water, can hardly be used to interpret the time delay of lowfrequency anomalies, especially in the case of a long time delay. This study investigates the generation mechanism of the LFSs based on a field case in the Ordos Basin, which is situated in the North China block (Wang et al. 2005). The paper is organized as follows. Section 2 describes geological background, petrophysical characteristics and seismic LFSs of the study area. Some parts of the target horizon in the gas field show LFSs, but some parts do not have any LFSs. Section 3 details the related petrophysical model. Section 4 analyses the effects of the pore-fluid property, and the thickness and permeability of reservoirs on the frequency-dependent reflection coefficients that are related to the low-frequency shadows.
Finally, section 5 discusses mechanisms of strongamplitude low-frequency shadows. The pore-fluid properties may change the fluid-pressure diffusion, which in turn affects the fluid-pressure equilibration between the adjacent layers, resulting in different frequency-dependent attenuation phenomena (White 1975). Gas-bearing reservoirs may cause strong absorption. When seismic wavelengths are much greater than the thickness of the reservoir, Biot's slow P-wave will occur and there will be equilibration of the fluid pressure between layers. The small thickness of reservoirs implies interference effects, which affect the frequency-dependent reflection coefficient and the arrival time of the seismic wave.

Geological background
The study area is located on the north Yi-Shaan slope in the Ordos Basin, China (figure 1). The typical tectonic form is a huge but gentle monocline dipping SW with a 1-2 o slope angle (Yang et al. 2012). There is a large-scale continuous lithological gas reservoir, the Sulige gas field (Yang et al. 2017). This gas field was discovered in year 2000 when the first well yielded a gas flow of 120 × 10 4 m 3 per day from Member H 8 of the Permian Lower Shihezi (LS) formation (Fu et al. 2019). The pay zone of this gas field includes Member H 8 of LS and Member S 1 of the Permian SX formation, beneath LS ( figure 1c).
The regional geology shows that Member H 8 of LS mainly developed fluvial-delta facies, which include underwater distributary channels and braided channel microfacies (Fu et al. 2008). The gentle tectonic setting results in the large-scale continuous distribution of super-imposed shaly sand reservoirs. The single high-quality sand body is small in scale and is stacked with other sand bodies to form a large, and effectively complex, sand body.
The lithology of Member H 8 of LS includes quartz sandstone and lithic quartz sandstones. Their mineral components are primarily quartz, feldspar, calcite and some clay minerals. The high porosity sand body is a favourable place for natural gas storage. It is believed that the gas in reservoirs comes from the coal seam strata beneath Member H 8 of LS (Fu et al. 2019). The water-bearing sand bodies are widely distributed, while the gas-bearing sand bodies have a limited distribution in the Sulige area. Figure 2 displays a seismic profile from west to east in the centre of the gas field. Overall, the structure of Member H 8 of LS is a westward-tilting monocline. The seismic event at the bottom of Member H 8 (labelled H8-b) is characterized by a high amplitude, an intermediate frequency and an intermediate continuity. The most significant and strongest event (labelled B1) in figure 2 is a reflection from the coal seam strata, which is the gas source for the reservoirs in Member H 8 of LS.
There are four vertical wells drilled through the seismic profile shown in figure 2, among which wells A, C and D are gas-producing wells, and well B is a water-producing well. It is noteworthy that the reflected energy from the gas-saturated reservoir (indicated by the green solid circle) within Member H 8 is slightly stronger than that of the watersaturated reservoir (indicated by the blue solid circle). Member H 8 at well C is structurally higher compared to the other wells, whereas Member H 8 at well A is structurally lower. Well B is located in the water-saturated part of Member H 8 that is deeper than the gas-saturated part of the formation at well C and well D. Obviously, structure is not an important factor when predicting the gas-bearing reservoirs in this area. The existence of gas reservoirs is controlled not only by the distance between the gas source rocks and the reservoirs (Meng et al. 2016), but also by the differences in the physical properties of the reservoir at different fluid saturations. location on the north Yi-Shaan slope. The typical tectonic form is a huge but gentle monocline dipping SW with a 1-2 o slope angle. (c) Composite stratigraphic column from a typical well. Corresponding to the regional geology, BX is the Ben-xi formation, TY is the Tai-yuan formation, SX is the Shang-xi formation, LS refers to the Lower Shihezi formation and US is the Upper Shihezi formation.

Petrophysical characteristics
Borehole data suggest that the sand reservoirs within H 8 Member of LS have a discrete vertical distribution of sands and shales, as outlined in figure 3. Vertically, the gas-and water-bearing sand zones are poorly connected and they are distributed in isolation. The water-gas interface in Member H 8 of LS is not uniform, and even within a single sand reservoir, because it is difficult to differentiate between gas and water (Meng et al. 2016). Figure 4a presents a core segment of Member H 8 of LS taken from well C. There are several horizontal thin layers with thicknesses in the order of a few centimetres or more. These thin layers may have been saturated with different fluids at the same time due to changes in the lithologies, water saturation and pore structures of the surrounding rocks. The reservoir within Member H 8 is characterized by complex pore throat structures, low porosity and permeability, and high saturation. The formation water in these wells occurs as free water, irreducible water and retained water and coexists with the gas. The complexity of the gas-water distribution within Member H 8 makes seismic identification of the gas reservoirs very difficult.
As indicated in figure 4b, the pore structure of the reservoir is dominated by secondary pores, such as dissolution pores. The results of fluid displacement experiments on rock samples from Member H 8 of LS in the study area show that the relationship between the P-wave velocity and the water saturation can be interpreted using Brie's empirical relation (Brie et al. 1995) with exponent of 1-3. This is very close to the Voigt Average approximation, indicating that the fluid distribution within the rocks is mainly characterized by patchy saturation. According to the statistics of the logging data for Member H 8 of LS at these wells, the difference between the impedance of the reservoirs and surrounding rocks is very small, which makes it very difficult to predict the locations of the reservoirs using conventional seismic inversions. It is well known that fluid saturated in the form of patchy saturation is one of the main causes of seismic wave-induced fluid flow, which leads to frequency-dependent seismic attribute (such as seismic attenuation) variation with pore-fluid changes in the seismic frequency band.

Low-frequency shadows associated with the gas reservoir
To investigate the seismic response of the fluid-saturated reservoir within Member H 8 of LS at different frequencies, we generated time-frequency spectra by applying the S-transform to the seismic data (figure 2). The time-frequency spectra are generated by decomposing the seismic waveform that is centred at each time sample into its constituent frequency. The spectral variation pattern at each time sample can be used for reservoir characterization (Wang 2012). We analysed the seismic response of the fluid-saturated reservoirs spectrally at each frequency with a relatively high temporal resolution. Figure 5 displays the spectra at 8, 12 and 45 Hz. At the frequency of 8 Hz, we observed a series of obvious amplitude anomalies beneath the gas-bearing reservoirs of Member H 8 at wells A, C and D. However, there are no visible amplitudes at well B. The low-frequency amplitudes increased at 12 Hz. These amplitude anomalies correspond to the gas-saturation zone of the reservoirs and can be interpreted as the lowfrequency shadows. The reason for this amplitude behaviour is partly due to the existence of natural gas in the pore spaces and the thinly layered configuration. The reflected energy of the gas-bearing reservoirs is less attenuated at low frequencies than at high frequencies, and the thin layers may delay the seismic response of the gas-bearing reservoirs. As expected, the anomalies related to the LFS are reduced at a frequency of 45 Hz. The amplitudes become weaker, as shown in figure 5c, and the anomalies only represent stratum interfaces.

The poroelastic model
To interpret LFSs related to gas reservoirs, we constructed a reservoir model based on the well logs and the core data. We have conventional wireline logs for the four boreholes, including gamma ray (GR), spontaneous potential (SP), bulk density (RHOB), sonic (DT), compensated neutron logging (CNL) and deep, medium and shallow electrical resistivity logs. In addition, the acoustic data from the multipole array logging (MPAL) are available for well C.
Thirty-five core samples are available for Member H 8 of the Permian LS formation in the study area and all samples are shaly sandstones. Core measurements were conducted to determine the porosity and permeability values, and the P-and S-wave velocities of the drained and undrained samples (Fang et al. 2019). The mineral contents of the samples were determined using the results of X-ray diffraction (XRD), and the pore structures were characterized using scanning electron microscopy (SEM).
Before the well data were interpreted, the GR, RHOB, CNL, DT and deep electrical resistivity curves of each well were normalized to help maintain good consistency among these wells, resulting in accurate estimations of the rock properties from the well data. Nine model parameters were needed for the poroelastic simulations, i.e. three poroelastic constants, three petrophysical parameters and three fluid property-related parameters (Li & Rao 2020). The initial step in constructing the reservoir model in our study was to generate a lithological volumetric model of the target zone. A five-component model was used, i.e. quartz, calcite, feldspar, clay and porosity, based on the information that was available from the XRD results and the logging data for the wells. The GR, RHOB, CNL, deep electrical resistivity and SP curves were used to construct the component model. An optimization algorithm was used to obtain the mineralogical volumetric fractions and the porosity from the logging data. The results of the multi-mineral decomposition were constrained and calibrated using the mineral composition that was determined by the results of the XRD. Figure 6 shows the constructed reservoir model based on the logging data for well C and the laboratory-derived measurements. The porosity of the model (track 5 in figure 6) was determined using the results of the optimization solution and was calibrated using the porosity that was derived from the plot of the RHOB versus the CNL porosities. In addition, a Kozeny-type equation was used to develop the empirical relationship between the permeability and porosity using the measured data, and we obtain where is the permeability in Darcy, and is the porosity (decimal). This relationship was used to calculate the permeability from the porosity curves. The solid bulk modulus of the model was computed using Hill's averaging based on the elastic properties and the volumetric fractions of the mineral compositions that were determined using the optimization solution, as shown in track 1. The solid bulk density was calculated using the arithmetic average of the density and volumetric fractions of each composition. The drained bulk and shear moduli were estimated using the equations given by Pride et al. (2002): where K d and K s are the bulk moduli of drained and grain materials, respectively, d and s are the shear moduli of the drained and solid grain materials, respectively, and a is a free parameter that depends on the pore structure, the rock consolidation and the ratio K s ∕ s . We estimated parameter a using equations (2) and (3) and the porosity data (track 5), RHOB curve and S-velocity of the array acoustic logging data. The estimated K d and d are shown in tracks 2 and 3 in figure 6. The well logging data were sampled at 0.125-m intervals. To determine the heterogeneity of the lithologies in the target zone, we applied a proximal support vector machine classifier algorithm (Mangasarian & Wild 2006) to the well data. Then, the classification results were calibrated using the lithologic information obtained from the logging data of the borehole drill cuttings.
In the Sulige gas field, the pore pressure of Member H 8 was in the range of 25-31 MPa, and the temperature of the target zone was about 102°C. Based on the results of the formation water analysis, the salinity is about 55 000 mg l −1 . Then, we estimated the fluid bulk moduli of the formation water and gas to be 2.25 and 0.02 GPa, respectively, the fluid densities to be 1000 and 150 kg m −3 , and the fluid viscosities to be 0.001 and 0.00002 Pa s −1 , respectively. Finally, the values of fluid bulk moduli, densities and viscosities at each sampled depth were assigned based on these parameters and the lithofacies results, and the water saturation derived from the electrical resistivity data, as shown in tracks 7-9.

The effect of seismic frequency, reservoir thickness and pore-fluid property
To understand the LFS phenomenon, we investigated the effect of the frequency, the reservoir thickness and the porefluid properties on the seismic response. We numerically modelled the low-frequency seismic response of the target zone using the poroelastic simulation method (Li & Rao 2020). We calculated the seismic responses of the established model, which is saturated with different fluids. The seismic response is a convolution of the total P-wave reflection coefficients convoluted with a normally incident Ricker wavelet, which has different peak frequencies (Wang 2014(Wang , 2015. Figure 7a shows the low-frequency (10 Hz) synthetic seismogram for the reservoir model (figure 6), which is composed of layers of 0.125 m in thickness. As displayed in figure 6, the total thickness of the reservoir model is about 31 m, which is about one-eighth of the wavelength. Thus, the seismogram in figure 7a depicts the integrated seismic response of the fluid-saturated thin layers in the reservoir model. There is also a delayed seismic response immediately after the primary P-wave reflections (the red ellipse). This delayed seismic response is likely caused by the mode conversion of the fast-to-slow P-waves. The delayed time of the seismic responses is comparable with the delay of LFS in well C (figure 5b).
For comparison, figure 7b is the synthetic seismogram that was generated by using the same parameters, but with the pore space of the reservoir model fully saturated with water, instead of the fluid in the field case. We can see that the seismic response from the hypothetical water-saturated model has a shorter duration and a fast energy attenuation at the tail of the seismic response. As expected, the amplitude of the seismogram for the water-saturated model is smaller than that for the gas-saturated model in the field case.
The numerical results (figures 7a and 7b) indicate that the seismic responses of thinly layered fluid-saturated reservoirs can produce a significant gas shadow at low frequencies. This LFS can be attributed to the existence of converted fast-to-slow P-waves. To illustrate the reasonableness of this conjecture, we carried out a simulation (figure 7c) using the same model parameters but having a Ricker wavelet with a higher peak-frequency (45 Hz). The results are shown in figure 7c. Now, if we compare the low-frequency synthetic trace of figure 7a to the high-frequency synthetic trace of figure 7c, we can see clearly that the low-frequency seismic responses in figure 7a are significantly prolonged, and the delayed low-frequency amplitudes are significantly enhanced. This is expected, since the energy of the converted fast-toslow P-wave propagating within the fluid-saturated media is exponentially attenuated, and the converted slow P-wave in the low-frequency range is more developed than that in the high-frequency range.
We regard the difference in the duration of the seismic response between the low-frequency seismogram (figure 7a) and that shown in figure 7c as the delay time. The delay of the low-frequency seismic response is about 30 ms, which is close to the delay for well C ( figure 5b).
Moreover, the reservoir thickness plays an important role in the total P-wave reflections. To investigate the effect of reservoir thickness on the delay of LFSs, we defined the model parameters with a uniform layer thickness of 3 m. The average values of the physical parameters of each layer were determined using the following equations: where the superscript j is the layer number, and ⟨ ⋅ ⟩ is the averaging operator across the calculated thickness. The poroelastic method was applied to the redefined reservoir model using a Ricker wavelet of 10 Hz. The result shown in figure 7d reveals that the converted slow P-wave effects died out as the layer thickness is increased. Thus, it only has a primary seismic response of the fluid-saturated model without the presence of sidelobes, caused by the interference of the converted slow P-wave. Figure 7d shows a strong contrast to the synthetic trace shown in figure 7a. Reflection coefficients of the total converted fast-to-slow P-waves. Blue and green curves indicate the saturation with water and gas, respectively.

Discussion
The observed seismic amplitude difference between gas-and water-saturated models in the previous section is consistent with the numerical results (figure 8) for the layered model saturated with different fluids. The layered reservoir model consists of a reservoir sandwiched between two identical background layers. The physical properties of these layers (listed in Table 1) correspond to the typical rock and fluid parameters of the seal and reservoir within Member H 8 in the study area. We set the pore space of the reservoir to be fully saturated with water and gas, respectively, and the thickness of the reservoir to be 1 m. The fluid properties are listed in Table 2. The calculation method is presented in the paper of Li & Rao (2020). The calculation results (figure 8) indicate that both the total P-wave reflection coefficients (R t p ) and the total converted slow P-wave reflection coefficient (R t cp ) increase with increasing frequency. R t p decreases as the fluid bulk moduli increase, but R t cp increases as the fluid bulk moduli increase. The absolute values of R t cp is three orders of magnitude smaller than that of the R t p for the model saturated with different fluids.
The layer thickness of reservoirs also has an effect on the seismic response (Yang et al. 2016). We set the pore space to be saturated with gas, and the thicknesses of the reservoir to be 0.3, 3.0 and 15 m, respectively. The calculation results (figure 9) indicate that layer thickness significantly influences both R t p and R t cp . For the model with a greater thickness of reservoir, the influence of that thickness is very weak. The converted slow P-wave effect is probably negligible for models in which the thickness of the reservoir is greater than 3.0 m, due to the high attenuation of the slow P-wave. However, if the thickness of the reservoir is thin enough, such that the converted slow P-wave generated at the top of the reservoir can reach the bottom and reflect back to the top of the reservoir again, then the converted slow P-wave may significantly influence the total reflection coefficient. There is a peak in R t p (figure 9a) corresponding to the model with a thickness of 0.3 m, which may be explained by the additive effect arising from the positive interference of the slow P-wave.
Finally, a high heterogeneous factor that affects the behaviour of slow P-wave is the rock permeability. The permeability of sandstone reservoirs can span seven orders of magnitude. Increasing the permeability often results in the amplitude of a slow P-wave increase (Bouzidi & Schmitt 2009). We considered the effect of the permeability of the reservoir, which is fully saturated with gas. Figure 10 displays R t p and R t cp when the permeability of the reservoir is 50, 500 and 1500 mD. The permeability of the reservoir 9 Figure 9. Effect of reservoir thickness on reflection coefficients as a function of frequency. (a) Reflection coefficients of the total P-wave. (b) Total converted fast-to-slow P-waves. The reservoir thickness is set to 0.3 m (red curves), 3 m (green curves) and 15 m (blue curves). Figure 10. Effect of reservoir permeabilities on reflection coefficients, as a function of frequency. (a) Reflection coefficients of the total P-wave as a function of the frequency. (b) The total converted fast-to-slow P-waves as a function of the frequency. The permeability of the reservoir is 50 mD (blue curves), 500 mD (red curves) and 1500 mD (green curves). affects both R t p and R t cp because it changes the behaviour of the converted slow P-wave. At higher frequencies, increasing the permeability of the reservoir increases both R t p and R t cp , since the amplitude of the converted slow Pwave generated at an interface increases with increasing permeability.
Numerical results show that the R t p of porous layered media is frequency dependent and is affected by the pore-fluid properties, as well as the thickness and rock properties of the layer. R t p increases when the permeability or frequency increases. The converted slow P-wave may change the tuning effect of the R t p of thinly layered media, especially for a model with a layer thickness on the order of tens of centimetres. Although R t cp is 3-5 orders of magnitude smaller than the R t p in our examples, the existence of the converted slow P-wave significantly influences the low-frequency P-wave response of thinly layered media.
In the previous section, we have shown the seismic responses of the same reservoir model at different frequencies (figure 7). To highlight the late low-frequency anomalies that are used for seismic fluid identification, it is convenient to use their differences. An effective way to reduce the uncertainty in the seismic predictions is to use the spectral attributes at different frequencies, since the delay time mainly depends on the existence of the converted slow P-wave. Hence, we calculated the difference of spectral amplitudes between 12 Hz (figure 7b) and 45 Hz (figure 7c). The result is shown in figure 11. 10 Figure 11. Amplitude-difference attribute, obtained from spectral subtraction of 12 and 45 Hz iso-frequency data. The high values of the attribute (warm colours) fit gas-producing wells A, C and D. The zone on the left of Well C beneath Member H 8 shows an anomaly, which might be associated with the gas-bearing reservoir.
As demonstrated in figure 11, seismic data at the gasproducing wells A, C and D generate obvious LFSs. The energy of the low-frequency anomalies has been enhanced, and the distribution of the low-frequency energy has been more concentrated. Nevertheless, water-producing well B does not exhibit any low-frequency anomalies. Figure 11 also clearly shows that there is a high attribute zone immediately beneath Member H 8 to the left of well C. Based on the relationship between the existence of LFSs and the gas-bearing reservoirs, it is likely that this low-frequency anomaly zone is associated with the gas reservoirs (as marked in a red ellipse).
The energy and delay time of the LFS of each gasproducing well are different (figure 11). Our numerical results cannot quantitatively explain the exact characteristics of the LFS, but they can be used to qualitatively interpret the relationship between the LFS and the gas-bearing reservoirs. Although significant efforts have been made to build a real reservoir model and to deliberately simulate the lowfrequency energy anomalies of the model, it is difficult to establish a standard LFS pattern for the study area. As previously demonstrated, the effect of the rock properties, such as the permeability and elastic contrasts, affects the converted slow P-wave and, in turn, the LFS. As every rock has its own history, different gas-bearing reservoirs may produce different patterns of low-frequency energy anomalies.

Conclusions
We have simulated low-frequency responses of a real gasbearing reservoir, using a poroelastic simulation method, to quantify the characteristics of seismic low-frequency shadows. Seismic data, well logging and core measurements in the study area indicate that the gas-bearing reservoir of the study area is highly heterogeneous in both the horizontal and vertical directions. The water-gas interface in the reservoir within the formation is not uniform. The existence of pore-fluid patchy saturation within the reservoir forms the condition for frequency-dependent seismic response analysis. The simulated result has been in good agreement with seismic data and confirms the validity of these LFSs in the field case. The investigation of the factors affecting the characteristics of LFS has shown that both pore-fluid properties and reservoir thickness are the key elements for interpreting the LFSs, and the delay time mainly depends on the existence of the converted fastto-slow P-waves. In practice, however, a detailed structure of the rock elastic properties, including the fluid content and the porosity and permeability at the scale of centimetres, needs to be defined, for quantitatively predicting the amplitude intensity and delay time of the low-frequency shadows.