Characterization of a Volcanic Gas Reservoir Using Seismic Dispersion and Fluid Mobility Attributes

Seismic dispersion and ﬂ uid mobility attributes are used to characterize a volcanic gas reservoir in the Songliao Basin of China. A rock physics model is constructed to describe poroelastic behaviors associated with heterogeneous ﬂ uids saturation within the volcanic gas reservoirs, where velocity dispersion and attenuation of propagating waves are attributed to the wave-induced ﬂ uid ﬂ ow described by the patchy saturation theory. Modeling results indicate that the frequency-dependent bulk modulus at the seismic frequency is more sensitive to gas saturation than the P -wave velocity dispersion. Accordingly, a new inversion method is developed to compute bulk-modulus-related dispersion attribute D K for improved characterization of volcanic gas reservoirs. Synthetic tests indicate that D K is more sensitive than traditional P -wave dispersion attribute D P to the variations of reservoir properties. The high value of dispersion attribute D K indicates the volcanic gas reservoirs with high porosity and gas saturation. At the same time, ﬂ uid mobility attribute FM can discriminate the volcanic gas reservoir as D K . Field data applications illustrate that D K and FM exhibit anomalies to the gas zones in the volcanic gas reservoir on the cross-well section. However, D K is more robust than FM to identify favorable zones on horizontal slices for speci ﬁ c target layers. Overall, rock physical modeling provides insights into the poroelastic behaviors of volcanic gas reservoirs, and inversion for seismic dispersion attribute D K improves hydrocarbon detection in the volcanic gas reservoir.


Introduction
As a type of essential unconventional reservoir, the volcanic reservoir plays a significant role in hydrocarbon production. Amongst, the volcanic gas reservoir of the Songliao Basin in northeast China contributed an increasing amount of gas production in recent years. This volcanic gas reservoir is located in the formations related to multiple-period volcanic eruptions, representing a heterogeneous stratigraphic structure with very low porosity and permeability. Complex microstructure makes it challenging to get insights into rock physical properties and associated seismic responses of the volcanic gas reservoir. Meanwhile, it is demanding to develop practical seismic techniques to discriminate favorable volcanic gas reservoirs based on rock physical modeling.
Seismic techniques such as amplitude variations versus offset (AVO) are usually used to identify lithology and fluids in a hydrocarbon reservoir. Various linear approximations of reflection coefficients were proposed based on the Zoeppritz equation for AVO inversion. Aki and Richards [1] derived an approximate formula of the PP -reflection coefficient in terms of P-and S-wave velocities and density from the Zoeppritz equation. Shuey [2] further proposed an approximate formula related to Poisson's ratio. Smith and Gidlow [3] developed an approximate expression represented by the P-and S-wave velocities. Then, Gray [4] gave a linearized approximation using the bulk and shear moduli. Furthermore, Russell et al. [5] proposed a generalized AVO representation to calculate fluid properties by incorporating the poroelastic theory.
With the progress in rock physical methods, velocity dispersion and attenuation associated with wave-induced fluids flow in a porous medium with complex microstructures were investigated [6][7][8][9][10][11][12][13]. Chapman et al. [14] found that the reflection coefficient for hydrocarbon reservoirs represents strong frequency-dependence and thus suggested characterizing reservoirs using seismic dispersion anomaly.
Accordingly, Wilson et al. [15] proposed the frequencydependent AVO (FD-AVO) scheme using the linearized approximation of Smith and Gidlow [3] for improved identification of gas reservoirs. Zhang et al. [16] proposed the FD-AVO inversion scheme based on Shuey's approximation. Cheng et al. [17] proposed identifying tight gas sandstone reservoirs using frequency-dependent attributes regarding intercept, gradient, fluid factor, and pseudo-Poisson's ratio. Sun et al. [18] investigated the sensitivity of the P-wave velocity dispersion attribute to fluid types in hydrocarbon prediction for a carbonate reservoir. Han et al. [19] combined the azimuth AVO inversion with spectral decomposition and used anisotropic P-wave dispersion attribute to discriminate fractures and fluids. Du et al. [20] tested the influence of frequency range on FD-AVO inversion. Chen et al. [21,22] represented a case study of the FD-AVO inversion for a sandstone hydrocarbon reservoir and proposed a dispersion attribute inversion using poststack seismic data. Pang et al. [23] and Liu et al. [24] applied the FD-AVO method to identify sweet spots in shale gas reservoirs. At the same time, some researchers have made efforts to involve sophisticated spectrum decomposition techniques in the FD-AVO inversion. For instance, the variational modal decomposition [25] and the inversion spectrum decomposition [26] were applied as the spectrum decomposition methods to improve inversion resolution.
Besides dispersion attributes obtained from FD-AVO inversion, many other frequency-related seismic attributes have also been proposed for identifying hydrocarbon resources, such as high-frequency attenuation [27], lowfrequency shadows [28], and fluid mobility [29,30]. The fluid mobility (FM) attribute was proposed based on a simplified approximation of reflection coefficients at the low frequency in poroelastic media. Cai [31] systematically studied the relationships between fluid mobility and reservoir permeability and hydrocarbon production. Chen [32] proved that fluid mobility could be positively correlated with the derivative of the reflection coefficient for the frequency and suggested using the instantaneous spectrum of seismic data in applications. At the same time, Wang et al. [33] and Luo et al. [34] predicted gas reservoirs using the fluid mobility derived from the low-frequency component of the marine seismic data.
At present, hydrocarbon identification using seismic dispersion attributes for volcanic reservoirs has not been investigated sufficiently. Moreover, the rock physical mechanism of velocity dispersion in volcanic rocks remains not so clear. In this study, we propose a new method to compute seismic dispersion attributes for improved identification of the volcanic gas reservoir. We build a rock physics model to investigate the poroelastic properties of volcanic rocks related to wave-induced fluid flow due to heterogeneous fluid distribution. Next, seismic dispersion attribute inversion methods are tested using synthetic data for theoretical models with different reservoir properties. Finally, the proposed methods are applied to field data to characterize a volcanic gas reservoir in the Songliao Basin, northeast of China.

Bulk-Modulus-Related Seismic Dispersion Attribute.
Gray [4] proposed the approximation of the reflection coefficient in terms of the bulk modulus and the shear modulus based on Aki and Richards [1] where K and μ are the bulk and shear modulus and ρ denotes density; the symbol Δ represents the difference of corresponding properties across an interface; θ is the incidence angle; γ sat = ðV P /V S Þ sat indicates the ratio of the P-and S-wave velocity of the saturated rock. We simplify equation (1) to by setting AðθÞ = ð1/4 − 1/3γ 2 sat Þ sec 2 θ, BðθÞ = 1/γ 2 sat ð1/3 sec 2 θ − 2 sin 2 θÞ, and CðθÞ = 1/2 − 1/4 sec 2 θ.
In equation (2), we assume that density ρ is independent of frequency and then obtain the frequency-dependent AVO representation where the bulk K and shear modulus μ vary with the seismic frequency. Then, we expand equation (3) on the frequency using Taylor's first-order expansion where f 0 is the reference frequency within the seismic frequency. By defining the reflection coefficient at the reference frequency we write equation (4) as In equation (6), we define the modulus-related dispersion attributes D K and D μ , which have the form of Then, we rearrange equation (6) as where Equation (10) can be expressed in the following form of the matrix for field data applications The matrix involves a prestack seismic angle gather with n incidence angles and m frequencies obtained from spectrum decomposition at each time interval.
In equation (11), we incorporate wavelet spectrum ½Wðf 1 Þ, Wðf 2 Þ, ⋯, Wðf m Þ to consider the overprinting effect of the incidence wavelet on the decomposed spectra of the prestack angle gathers ΔSðθ, f Þ. For synthetic tests, the wavelet spectrum is calculated from the incidence wavelet in the forward modeling. For field data applications, it can be computed from the extracted wavelet.
We simplify the equation (11) in terms of where d represents the decomposed spectrum ΔSðθ, f Þ; G indicates the coefficient matrix; m is the vector of D K and D μ to be computed. In equation (12), m can be solved by the damped leastsquares inversion method, which has the form of where ε is the damping factor; I is the identity matrix.

Fluid Mobility Attribute.
Silin et al. [30] researched seismic reflection from a boundary between elastic and fluidsaturated porous media by incorporating poroelastic theory.
The frequency-dependent reflection coefficient for a normal incidence P-wave was represented as where κ is permeability and η is fluid viscosity; i is the imaginary unit; ρ b is rock bulk density; R 0 and R 1 are real and frequency-independent terms related to porosity and elastic properties of fluids and solid matrix. Chen et al. [32] and Wang et al. [33] used the factor FM = κ/η as an indicator of fluid mobility and then expressed equation (14) as By taking derivative to the frequency for equation (15) as follows and setting C = ffiffiffiffiffiffiffi π/2 p R 1 ð1 + iÞ ffiffiffiffi ffi ρ b p , the fluid mobility attribute can be expressed as 3 Lithosphere Equation (17) shows that FM is proportional to the squared derivative of the reflection coefficient to the frequency at a specific frequency. Chen et al. [32] and Wang et al. [33] suggested representing the magnitude of the reflection coefficient at a particular frequency f using decomposed spectrum Sð f Þ. Thus, fluid mobility attribute FM can be estimated by where Δf represents the frequency interval.

Rock Physics Model of Volcanic Gas Reservoir.
According to the geological analyses of the study area, the volcanic gas reservoir represents low porosity and permeability. Thus, the traditional Gassmann fluid-substitution theory does not work well for tight volcanic rocks. Meanwhile, the complex pore space of the volcanic gas reservoir comprises intergranular pores and induced microcracks as shown in Figure 1(a), which should be considered in rock physical modeling. In Figure 1(b), we built a rock physics model to calculate the properties of volcanic rock associated with complex microstructure and fluid flow. In the modeling workflow, we apply the Hashin-Shtrikman bounds (HSB) proposed by Hashin and Strikman [35] to calculate elastic properties of the solid matrix, which mainly consists of quartz and clay minerals. Then, we employ the self-consistent approximation (SCA) method given by Berryman [36] to put pores and microcracks into the solid matrix and obtain the elastic properties of the dry skeleton. Finally, we use the patchysaturated model given by White [37] to calculate poroelastic properties associated with heterogeneous fluid saturation in volcanic rock. Table 1 shows the properties of constituents for rock physical modeling. According to laboratory measurements on the volcanic rock, we assume volumetric fractions of quartz 80% and clay 20%, and the aspect ratios of pores and microcracks in the dual-pore model are 1 and 0.01.
The proportions of pores and microcracks in the total pore space are 0.89 and 0.11, respectively.
In the workflow in Figure 1(b), the HSB theory estimates the upper and lower elastic bounds K HS± and μ HS± where K and μ represent the bulk and shear modulus; f is the volumetric fraction of compositions. Subscripts 1 and 2 denote the two components that form solid matrix. Elastic moduli K s and μ s of the solid matrix are estimated by averaging the upper and lower bounds in equations (19) and (20).   Berryman [36] proposed the SCA theory to compute the effective bulk and shear moduli K * SC and μ * SC of a multiphase medium by where index i denotes the i-th constituent; x i represents volumetric fraction; K i and μ i are bulk and shear modulus of the i-th constituent; K * i and Q * i are the geometric factors of the inclusions put into the background material.
According to the patchy saturation theory of White [37], heterogeneous fluid saturation and wave-induced fluid flow are regarded as the mechanism to cause dispersion and attenuation of seismic waves propagating in poroelastic where V * P = ffiffiffiffiffiffiffi E/ρ p is complex P-wave velocity; E = K * SC + 4μ * SC /3 denotes the P-wave modulus; density of the patchysaturated medium is calculated by ρ = ð1 − φÞρ s + φρ f ; ρ s 6 Lithosphere and ρ f are the density of solid matrix and fluids, and φ is porosity.
In Table 2, we design three models for the volcanic gas reservoir with different porosity and gas saturation. From model 1 to 3, porosity and gas saturation increase simultaneously. Thus, model 3 represents a high-quality volcanic gas reservoir compared to the others.
Then, we calculate the frequency-dependent P-wave velocity V P and bulk modulus K as shown in Figures 2(a) and 2(b), respectively. Between model 1 and model 3, the change rate of V P is about 9.5%, and that of K is about 27%. Therefore, the bulk modulus K is more sensitive to porosity and gas saturation changes than the P-wave velocity V P . Thus, in this study, we propose the bulk-modulusrelated dispersion attribute inversion method for improved identification of volcanic gas reservoirs.

Synthetic Tests.
Based on logging data interpretation of the study area, we design theoretical models simulating the volcanic gas reservoir of the study area, as illustrated by the P-wave impedance I P in Figure 3(a). Layer thickness of the volcanic gas reservoir is set to 40 m. After the depth to time convention, the target layer of the volcanic gas reservoir ranges from about 100 to 120 ms. Curves with different styles indicate three models of the volcanic gas reservoir in Table 2. Using results calculated for the three models in Figure 2(a), we use P-wave velocity V P at seismic dominant frequency 30 Hz is to compute I P as shown in Figure 3

Lithosphere
We simulate prestack AVO responses using the propagator matrix method [39]. For P-wave incidence, the reflection and transmission coefficient vector r = ½R PP , R PS , T PP , T PS T is computed by where i P is the P-wave incident vector; A 1 and A 2 are the propagator matrices of upper and lower half-spaces; B α = Q N α=1 B α = Tð0ÞT −1 ðh α Þ ðα = 1, 2, ⋯, NÞ is the propagator matrix of N layers between the two half-spaces; the matrix T is the function of the thickness h α for each layer.
The propagator matrix method (PMM) is implemented in the frequency domain. Thus, it can be conveniently linked to the poroelastic models describing velocity dispersion and attenuation in the frequency domain.
After computing the frequency-dependent reflection coefficient R PP using equation (25), reflection spectrum U f can be calculated for a given incidence wavelet with a spectrum W f Then, the reflected waveforms u t in the time domain can be obtained using the inverse Fourier transform  Figure 3(b). We calculate the seismic dispersion attribute D K and fluid mobility attribute FM, as shown in Figures 3(f) and 3(g),  9 Lithosphere respectively. We can see that the high D K and FM correspond to the volcanic gas reservoir with high porosity and gas saturation, as denoted by model 3. For comparison, we show the traditional P-wave dispersion attribute D P in Figure 3(h).
For quantitative comparison, we pick peak values of various attributes D K , D P , and FM for the target layer in Figures 3(f)-3(h) and display them in Table 3. For the magnitude of dispersion attributes between model 1 and model 3, the change rate is about 34% and 23% for D K and D P , respectively. Thus, D K is more sensitive than D P to the variations of reservoir properties. Meanwhile, the change rate of FM is about 83% between model 1 and model 3. It appears that FM is more sensitive compared to D K and D P . However, the robustness of reservoir property estimates using FM will be discussed later in this study.

Datasets.
Datasets used in this study are from the 3D seismic exploration for the volcanic gas reservoir of the Songliao Basin, northeast China. Figure 4 shows the logging data of well A in the study area, including gamma ray, Pwave velocity V P , S-wave velocity V S , density ρ, P-wave impedance I P , and porosity φ. Red lines indicate the strata formed by multiple-period volcanic eruptions. Periods 2 and 5, which have relatively high porosity and low density than adjacent layers, have been proven gas payzones.

Field Data Applications.
We then apply the methods to the field data. Figure 7 shows the sections of decomposed spectra at 5 Hz, 10 Hz, 20 Hz, 30 Hz, 40 Hz, and 50 Hz, calculated by implementing CWT on the poststack seismic section in Figure 5. The energy of seismic reflection concentrates around 10-30 Hz. Amplitudes from the top and bottom of the volcanic formations show substantial lateral variations on each frequency profile, implying possible changes in reservoir properties such as porosity and gas saturation. The frequency dependency of the seismic reflections revealed by the spectral decomposition is the basis to implement seismic dispersion attribute inversion and fluid mobility estimate.
We extract the wavelet in the time domain from the seismic dataset, as shown in Figure 8(a). Figure 8(b) displays the corresponding spectrum calculated using the FFT algorithm. The obtained spectrum in Figure 8(b) is incorporated in equation (11) to consider the wavelet effect in inversion.  Lithosphere Figure 9(b) illustrates the fluid mobility attribute FM computed using the method in Section 2.2. In equations (11) and (18), the involved frequency ranges from 5 to 50 Hz, covering dominant frequencies of seismic reflections. Sections in Figures 9(a) and 9(b) correspond to the poststack seismic profile shown in Figure 5. According to analyses using synthetic data in Section 3.2, high-value D K and FM indicate volcanic gas reservoirs with high porosity and gas saturation. Thus, the distributions of D K and FM in Figure 9 provide meaningful information for the characterization of the volcanic gas reservoir. Figure 10 shows the close-ups around well A for detailed analyses of D K and FM anomaly. Blue blocks indicate the locations of gas zones interpreted in well A. One gasproducing zone is just above the bottom of the volcanic formations, and the other is below and near the top boundary. We can see that both the D K and FM show high-value anomalies at the gas zones. However, it appears that compared to FM, D K represents more distinct and robust responses to the presence of gas zones.
Then, we apply the inversion methods to the 3D seismic dataset of the study area, with results shown in Figures 11  and 12. Figure 11 displays horizontal slices D K and FM for the top of the volcanic gas reservoir. The time window is 5-20 ms below the top horizon to cover possible locations of the gas zone. Figure 12 shows D K and FM for the bottom The upper gas zone shows distinct anomaly responses of seismic dispersion attribute D K , as indicated by the black dot in Figure 11(a). Anomaly responses of the lower gas zone are also apparent, as shown by the white dot in Figure 12(a). However, in Figures 11(b) and 12(b), FM does not estimate the upper and lower gas zones reasonably. Thus, D K plays a more robust indicator for identifying volcanic gas reservoirs than FM on the horizontal slices.
The reason for such different performances of D K and FM deserves further investigation. Possible explanations are when extending equation (2) to (3), the dispersion attribute D K inversion does not specify any velocity dispersion and attenuation mechanisms. Thus, D K exhibits more robustness than fluid mobility FM derived from Biot's fluid flow theory, which may be inapplicable to the volcanic rock in the study area.

Discussion and Conclusions
Estimating seismic dispersion attributes and fluid mobility is essential for better characterizing volcanic gas reservoirs. In this study, we built a rock physics model for the volcanic gas reservoir to describe poroelastic behaviors associated with heterogeneous fluid saturation in volcanic rocks. Velocity dispersion and attenuation of propagating waves are attributed to the wave-induced fluid flow described by the patchy saturation theory. The dual-pore model describes the complex pore structure of stiff round pores and soft microcracks in volcanic rocks. A modeling scheme was proposed by integrating poroelastic theory and the propagator matrix method (PMM) to simulate seismic signatures from the volcanic reservoir.
Rock physical modeling indicates that the bulk modulus and associated dispersion are more sensitive to gas saturation than the P-wave velocity. Thus, we developed an approach to calculate bulk-modulus-related dispersion attribute D K for improved reservoir characterization of the volcanic gas reservoir.
Synthetic tests indicate that the proposed dispersion attribute D K is more sensitive than the traditional P-wave dispersion attribute D P to the variations of reservoir properties. High-value D K corresponds to the volcanic gas reservoirs with high porosity and gas saturation. At the same time, fluid mobility attribute FM shows performance comparable to D K . The proposed method was applied to the volcanic gas reservoir of the Songliao Basin in northeast China. Field data applications illustrate that D K and FM represent apparent anomalies to the gas zones. However, D K is more robust than FM in identifying volcanic gas reservoirs on horizontal slices for specific target layers. The reason for the different performances of D K and FM in reservoir characterization deserves further investigation. Possible explanations are that the dispersion attribute D K inversion does not specify any poroelastic mechanisms. In comparison, fluid mobility FM is derived from Biot's fluid flow theory that may be inapplicable to the volcanic rock in the study area.
Overall, seismic rock physical modeling in this study provides insights into the poroelastic behaviors and corresponding seismic signatures of volcanic gas reservoirs. Inverted seismic dispersion attributes improves the hydrocarbon detection in the volcanic gas reservoir. In further studies, laboratory measurements on velocity dispersion and attenuation and more sophisticated rock physical modeling methods would help better understand the poroelastic properties of volcanic reservoirs. Based on this, more practical seismic attributes can be investigated for improved characterization of volcanic gas reservoirs.

Data Availability
No data were used, Available upon request or included within the article.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this article.