Seismic AVOA Inversion for Weak Anisotropy Parameters and Fracture Density in a Monoclinic Medium

In shale gas development, fracture density is an important lithologic parameter to properly characterize reservoir reconstruction, establish a fracturing scheme, and calculate porosity and permeability. The traditional methods usually assume that the fracture reservoir is one set of aligned vertical fractures, embedded in an isotropic background, and estimate some alternative parameters associated with fracture density. Thus, the low accuracy caused by this simplified model, and the intrinsic errors caused by the indirect substitution, affect the estimation of fracture density. In this paper, the fractured rock of monoclinic symmetry assumes two non-orthogonal vertical fracture sets, embedded in a transversely isotropic background. Firstly, assuming that the fracture radius, width, and orientation are known, a new form of P-wave reflection coefficient, in terms of weak anisotropy (WA) parameters and fracture density, was obtained by substituting the stiffness coefficients of vertical transverse isotropic (VTI) background, normal, and tangential fracture compliances. Then, a linear amplitude versus offset and azimuth (AVOA) inversion method, of WA parameters and fracture density, was constructed by using Bayesian theory. Tests on synthetic data showed that WA parameters, and fracture density, are stably estimated in the case of seismic data containing a moderate noise, which can provide a reliable tool in fracture prediction.


Introduction
In a variety of complex oil and gas reservoirs, fractured reservoirs widely exist in shale and igneous rock. Natural and induced fractures in reservoirs play an important role in determining fluid flow during production, and knowledge of the orientation and density of fractures is useful to optimize production from fractured reservoirs [1]. Areas of high fracture density may represent 'sweet spots' of high permeability, and it is important to be able to target such locations for infill drilling [2]. Recent developments in geophysics have revealed the viscous behavior of the underground strata, and demonstrated that the propagation of seismic waves has dispersion and attenuation, for fluid discrimination in porous rocks [3][4][5]. However, there are methods available to compensate for dispersion and attenuation. In the current implementation of the method, we also ignore the loss of seismic amplitude due to attenuation. Seismic anisotropy is defined as the dependence of seismic velocity upon angle. P-waves propagating parallel to fractures are faster than those propagating perpendicular to fractures [6]. For transversely isotropic (TI) media, empirical and analytical studies have shown that the presence of anisotropy can significantly distort conventional amplitude variation

Elastic Compliance Tensor of Fractured Medium
In an elastic medium that contains an arbitrary number of sets of fractures, with an arbitrary orientation distribution, by using the divergence theorem and Hooke's law, it can be shown that the elastic compliance tensor of the fractured medium can be written in the following form [29,30]: where S 0 is the compliance matrix of the medium (including the effects of pores, cracks, and stress, except for those fractures explicitly included in ∆S). The subscripts i, j, k, and l are the free index. The effective excess compliance ∆S ijkl , due to the presence of the fractures, with rotationally invariant shear compliance, can be written as [27]: Here, δ ij is the Krönecker delta, α ij is a second-rank tensor, and β ijkl is a fourth-rank tensor. If the two fracture sets are not orthogonal, but are perfectly aligned (within each set), the compliance tensors become [6]: cos ϕ 1 sin ϕ 1 0 cos ϕ 1 sin ϕ 1 cos 2 ϕ 1 0 0 0 0 sin 2 ϕ 2 cos ϕ 2 sin ϕ 2 0 cos ϕ 2 sin ϕ 2 cos 2 ϕ 2 0 0 0 0 sin 4 ϕ 1 sin 2 ϕ 1 cos 2 ϕ 1 0 0 0 2 sin 3 ϕ 1 cos ϕ 1 sin 2 ϕ 1 cos 2 ϕ 1 cos 4 ϕ 1 0 0 0 2 sin ϕ 1 cos 3 ϕ 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 sin 3 ϕ 1 cos ϕ 1 2 sin ϕ 1 cos 3 ϕ 1 0 0 0 4 sin 2 ϕ 1 cos 2 ϕ 1 sin 2 ϕ 2 cos 2 ϕ 2 0 0 0 2 sin 3 ϕ 2 cos ϕ 2 0 cos 4 ϕ 2 0 0 0 2 sin ϕ 2 cos 3 ϕ 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 sin 3 ϕ 2 cos ϕ 2 2 sin ϕ 2 cos 3 ϕ 2 0 0 0 4 sin 2 ϕ 2 cos 2 ϕ 2 where N V is the number of fractures per unit volume (fracture density), A is the average area of the fractures in the set, B T and B N are the (area-weighted) average specific tangential and normal compliances, and ϕ is the azimuthal angle between the fracture strike and the survey 1-axis. The subscripts 1 and 2 denote the 1st and 2nd fracture set, respectively. If the background medium is isotropic, or transversely isotropic, and there are at least two non-orthogonal vertical fracture sets, this leads to monoclinic symmetry of the fractured rock, as shown in Figure 1. are the (area-weighted) average specific tangential and normal compliances, and ϕ is the azimuthal angle between the fracture strike and the survey 1-axis. The subscripts 1 and 2 denote the 1st and 2nd fracture set, respectively. If the background medium is isotropic, or transversely isotropic, and there are at least two non-orthogonal vertical fracture sets, this leads to monoclinic symmetry of the fractured rock, as shown in Figure 1.  The additional compliance, due to two or more sets of aligned vertical fractures not aligned with the coordinate system, has the form [6]: When the fracture's additional compliance ∆S is small, one can directly calculate the stiffness matrix using the background stiffness matrix C 0 [2]: We assumed a polar anisotropic background. An isotropic background is a special case of this. Thomsen [31] proposed a simple elastic constant expression:

PP-Wave Reflection Coefficient and Generalized Anisotropy Parameters
In this study, the elastic contrast between the overburden and the reservoir was assumed to be small. In this situation, the plane-wave/P-wave reflection coefficient for a plane separating media, with arbitrary elastic symmetry and with WA, can be written in the form [32] 2 ∆ε x cos 4 φ + ∆ε y sin 4 φ + ∆δ z cos 2 φ sin 2 φ+2 ∆ε 16 cos 2 φ + ∆ε 26 sin 2 φ cos φ sin φ sin 2 θ tan 2 θ where R iso PP (θ) denotes the weak-contrast reflection coefficient at an interface separating two slightly different isotropic media, and the generalized Thomsen anisotropy parameters are given by where V p and V s are the P-and S-wave velocities of the background isotropic medium, respectively, and A ij = C ij ρ is the density-normalized elastic stiffness. Incorporating Equations (5) and (7) and substituting the density-normalized elastic stiffness of Equation (9), we obtained a new expression of the generalized Thomsen anisotropy parameters, for the monoclinic medium, in terms of WA parameters and the compliance tensors. These parameters are presented in Appendix A. Under the assumption of small compliance tensors and WA parameters, we neglected the product terms of the compliance tensors and the weak anisotropic parameters.
By substituting the generalized WA parameters in Equation (8) with Equation (A1), a P-wave reflection coefficient, in terms of WA parameters and the compliance tensors, was obtained: where the pseudo weak anisotropic parameter δ psd = M − µ + 2Mδ. The second term on the righthand side is only a function of offset because the background is VTI, which does not depend on the azimuth. Equations for sensitivities, F, are given in Appendix B.
In order to realize the direct inversion of fracture density, the linear P-wave reflection coefficient, in terms of weak anisotropic parameters and fracture density, could be obtained by combining Equations (3) and (4) and Equation (10).
where the sensitivities a(θ, φ) and b(θ, φ) for fracture density are given in Appendix C. The interpretation of the fracture reservoir showed that the fracture apertures within hydrocarbon reservoirs were 5 mm, and the majority were 1 mm or less [1,33]. The length of the open fracture is dependent on the formation pressure. The numerical analysis results showed that an open fracture with a width of 1 mm generally does not exceed 40 cm in reservoir depth. However, it is probably reasonable to think of the 40 cm value as an average distance between weld points. In that case, the fracture length can be 1 m long or 30 m long [33]. However, fracture orientations of shale can be easily obtained by field exposure or the formation microimager (FMI) log [19]. In this paper, under the assumption that the fracture width, radius, and strike were known, the WA parameters and fracture density could be directly estimated by Equation (11).

Bayesian Inversion for Weak Anisotropy Parameters and Fracture Density
Assuming that there are two sets of fractures, we produced a matrix of data equations for the case of m incidence angles and n azimuthal angles: From the above equation, it can be seen that elastic parameters of the isotropic background should be known for the accurate inversion of weak anisotropic parameters and fracture density. The V p , V s , and ρ values for all layers of interest were obtained from well logs using Backus averaging [34]. The forward problem has the simple form: d=WFm+n=Lm+n (13) where the forward operator L=WF, n is random noise added into clear data, and the model parameters are: and the wavelet matrix is: In Bayesian theory, the posterior probability density function (PDF) is calculated by the prior PDF and the likelihood function [35]: where p(m |d ) is the posterior PDF, p(m) is the prior PDF, and p(d |m ) is the likelihood function.
In the case of observed seismic data containing Gaussian random noise, the likelihood function p(d |m ) is given by where σ 2 e is the variance of the noise, N is the number of the input data samples, and G is the linear operator G = Lm.
Under the assumption of model parameters being independent of each other, the prior PDF p(m) was expressed as where C m is the variance of the model parameters.
According to Bayesian theory, we could convert the inversion problem by solving the minimum value of the objective function where R(m) = 1 2 m T C −1 m m, σ 2 d is the variance of data. By solving the derivative of the objective function with respect to the model parameters, we could obtain ∂m to zero gives the model parameters where super parameter µ h allowed us to estimate the model parameters along the tuning curve, mainly through experiments. In this paper, the derived P-wave reflection coefficient, after removing the background effect, was completely linear, thus allowing us to invert the model parameters from seismic data d.

Numerical Analysis
Synthetic seismic data were created by using interpretation results of well log data and vertical seismic profiling (VSP) data, which provided input for the AVOA inversion. It was assumed that the fracture medium had monoclinic symmetry, the azimuths of the two fracture sets were 10 • and 70 • , respectively, the average fracture width was 1 mm, the average fracture radius was 10 m, and the average area of the fracture was 0.02 m 2 . Figure 2 shows two sets of fracture densities from the well log interpretation results as model parameters for inversion. Figure 3 shows the well log and VSP interpretation results, among which the isotropic elastic parameters were used as the known parameters, to calculate the reflection coefficient of isotropic background and sensitivities.
where ( ) to zero gives the model parameters where super parameter h μ allowed us to estimate the model parameters along the tuning curve, mainly through experiments. In this paper, the derived P-wave reflection coefficient, after removing the background effect, was completely linear, thus allowing us to invert the model parameters from seismic data d.

Numerical Analysis
Synthetic seismic data were created by using interpretation results of well log data and vertical seismic profiling (VSP) data, which provided input for the AVOA inversion. It was assumed that the fracture medium had monoclinic symmetry, the azimuths of the two fracture sets were 10° and 70°, respectively, the average fracture width was 1 mm, the average fracture radius was 10 m, and the average area of the fracture was 0.02 m 2 . Figure 2 shows two sets of fracture densities from the well log interpretation results as model parameters for inversion. Figure 3 shows the well log and VSP interpretation results, among which the isotropic elastic parameters were used as the known parameters, to calculate the reflection coefficient of isotropic background and sensitivities.   In order to test the robustness of the proposed inversion algorithm, synthetic seismic data were created by the model parameters shown in Figures 2 and 3, as well as a 30 Hz Ricker wavelet, based on the derived PP-wave reflection coefficient, Equation (13), and the convolution model shown in Figure 4a. The incident angle had 10 angles between 0° and 45°, and the azimuths were 10°, 55°, 100°, and 145°.
Gaussian random noises (the signal-to-noise ratios (SNRs) were 5 and 2) were added into the synthetic data, as shown in Figures 4b,c. Comparisons between the inversion results (black) of the WA parameters, the fracture density, and the true values (red) are displayed in Figure 5. We can see that, in the case of the SNR being 5, the proposed inversion method could make a stable estimation of WA parameters and fracture density; however, when the SNR was 2, the inversion results of fracture density exhibited more instabilities. Random noise had some influence on the WA parameters and fracture density, but the overall agreement was good. In general, it can be concluded that the proposed inversion method is stable and valid.  In order to test the robustness of the proposed inversion algorithm, synthetic seismic data were created by the model parameters shown in Figures 2 and 3, as well as a 30 Hz Ricker wavelet, based on the derived PP-wave reflection coefficient, Equation (13), and the convolution model shown in Figure 4a. The incident angle had 10 angles between 0 • and 45 • , and the azimuths were 10 • , 55 • , 100 • , and 145 • .
Gaussian random noises (the signal-to-noise ratios (SNRs) were 5 and 2) were added into the synthetic data, as shown in Figure 4b,c. Comparisons between the inversion results (black) of the WA parameters, the fracture density, and the true values (red) are displayed in Figure 5. We can see that, in the case of the SNR being 5, the proposed inversion method could make a stable estimation of WA parameters and fracture density; however, when the SNR was 2, the inversion results of fracture density exhibited more instabilities. Random noise had some influence on the WA parameters and fracture density, but the overall agreement was good. In general, it can be concluded that the proposed inversion method is stable and valid.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 17 In order to test the robustness of the proposed inversion algorithm, synthetic seismic data were created by the model parameters shown in Figures 2 and 3, as well as a 30 Hz Ricker wavelet, based on the derived PP-wave reflection coefficient, Equation (13), and the convolution model shown in Figure 4a. The incident angle had 10 angles between 0° and 45°, and the azimuths were 10°, 55°, 100°, and 145°.
Gaussian random noises (the signal-to-noise ratios (SNRs) were 5 and 2) were added into the synthetic data, as shown in Figures 4b,c. Comparisons between the inversion results (black) of the WA parameters, the fracture density, and the true values (red) are displayed in Figure 5. We can see that, in the case of the SNR being 5, the proposed inversion method could make a stable estimation of WA parameters and fracture density; however, when the SNR was 2, the inversion results of fracture density exhibited more instabilities. Random noise had some influence on the WA parameters and fracture density, but the overall agreement was good. In general, it can be concluded that the proposed inversion method is stable and valid.     In the approach proposed, the fracture width, radius, and strike are required to be assumed. The dimensions of fractures are a priori functions of many factors in the medium and cannot be readily assumed in practice. In order to test the impact of this assumption (the dimensions assumed for the fractures) on the predicted results, the results of inversion, using the same synthetic data as used in Figure 4a, but with the forward operators, that fracture radius is 0.25 m and fracture width is 5 mm (wrong assumptions), are shown in Figure 6a and b, respectively. Comparing the inversion results with Figure 5a, we observed that, in the case of the fracture radius being 0.25 m, the incorrect assumption of fracture radius introduced larger errors into the prediction of WA parameters and fracture density. In general, the inverted model parameters captured the amplitudes and variability of the true models. When the fracture width was 5 mm, the incorrect assumption of fracture width almost had no impact on WA parameters; however, the accuracy and stability of the fracture density inversion was seriously influenced. More accurate information on fracture radius and width may help to improve the accuracy of WA parameters and fracture density inversion. In the approach proposed, the fracture width, radius, and strike are required to be assumed. The dimensions of fractures are a priori functions of many factors in the medium and cannot be readily assumed in practice. In order to test the impact of this assumption (the dimensions assumed for the fractures) on the predicted results, the results of inversion, using the same synthetic data as used in Figure 4a, but with the forward operators, that fracture radius is 0.25 m and fracture width is 5 mm (wrong assumptions), are shown in Figure 6a,b, respectively. Comparing the inversion results with Figure 5a, we observed that, in the case of the fracture radius being 0.25 m, the incorrect assumption of fracture radius introduced larger errors into the prediction of WA parameters and fracture density. In general, the inverted model parameters captured the amplitudes and variability of the true models. When the fracture width was 5 mm, the incorrect assumption of fracture width almost had no impact on WA parameters; however, the accuracy and stability of the fracture density inversion was seriously influenced. More accurate information on fracture radius and width may help to improve the accuracy of WA parameters and fracture density inversion. Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 17 (a) (b) Figure 6. Inversion results using synthetic data (Figure 4a), computed with the fracture radius at 10 m and fracture width at 1 mm. (a) Inversion results with the forward operator that the fracture radius was 0.25 m, (b) Inversion results with the forward operator that the fracture width was 5 mm.
In order to demonstrate the advantage of the proposed algorithm in fracture density estimation, we made use of the same synthetic data as used in the Figure 4a and estimated anisotropy through elliptical inversion of interval velocities. The resulting NMO velocity varied ellipsoidally with the azimuth. The anisotropy intensity was defined by the axial ratio. To compare with the true value of fracture density, the anisotropy intensity, multiplied by a factor, is shown in Figure 7. We can see that the accuracy of the proposed inversion, in relative terms, is greater than that of elliptical inversion. In order to demonstrate the advantage of the proposed algorithm in fracture density estimation, we made use of the same synthetic data as used in the Figure 4a and estimated anisotropy through elliptical inversion of interval velocities. The resulting NMO velocity varied ellipsoidally with the azimuth. The anisotropy intensity was defined by the axial ratio. To compare with the true value of fracture density, the anisotropy intensity, multiplied by a factor, is shown in Figure 7. We can see that the accuracy of the proposed inversion, in relative terms, is greater than that of elliptical inversion.

Conclusions
Assuming that the fractured reservoir has a monoclinic symmetry of VTI background, a new form of P-wave reflection coefficient, in terms of WA parameters and fracture density, was obtained by substituting the stiffness matrix coefficient of the VTI background, and normal and tangential fracture compliances. The new expression of the reflection coefficient for a monoclinic medium can avoid calculation errors due to the assumption of a simple model (such as an HTI medium). Additionally, the inversion parameters contain WA parameters, which avoid the assumption that the parameters of the VTI background are known in the traditional form of a P-wave reflection coefficient. Finally, we achieved the direct inversion of WA parameters and fracture density, and avoided the intrinsic errors introduced by indirect substitution. The synthetic test demonstrated that this method can be used to accurately estimate WA parameters and fracture density, which significantly benefits further calculation of reservoir porosity or permeability, the prediction of the "sweet spot", and the evaluation of a reservoir's reconstruction.
In the approach proposed, the fracture width, radius, and strike are required to be assumed. Incorrect assumptions of fracture width and radius will yield erroneous results, for WA parameters and fracture density, in practice. A priori knowledge about fracture width is important in the prediction of fracture density. However, the interpretation of fracture reservoirs shows that the fracture apertures in reservoir depth are about 1 mm, and have little variation. Funding: This research received no external funding.

Acknowledgments:
We appreciate the helpful comments and suggestions from editors and two anonymous reviewers that led to the substantial improvement of this article. We are grateful for the financial support of the following institutions: the Science and Technology Cooperation Project of the CNPC-SWPU Innovation Alliance and Open Projects Fund of the State Key Laboratory of Oil and Gas Reservoir Geology and Exploitation (PLN201733), and the National Natural Science Foundation of China (NSFC 41204101).

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
In terms of the specific compliance matrices, generalized anisotropy parameters are defined as below:

Conclusions
Assuming that the fractured reservoir has a monoclinic symmetry of VTI background, a new form of P-wave reflection coefficient, in terms of WA parameters and fracture density, was obtained by substituting the stiffness matrix coefficient of the VTI background, and normal and tangential fracture compliances. The new expression of the reflection coefficient for a monoclinic medium can avoid calculation errors due to the assumption of a simple model (such as an HTI medium). Additionally, the inversion parameters contain WA parameters, which avoid the assumption that the parameters of the VTI background are known in the traditional form of a P-wave reflection coefficient. Finally, we achieved the direct inversion of WA parameters and fracture density, and avoided the intrinsic errors introduced by indirect substitution. The synthetic test demonstrated that this method can be used to accurately estimate WA parameters and fracture density, which significantly benefits further calculation of reservoir porosity or permeability, the prediction of the "sweet spot", and the evaluation of a reservoir's reconstruction.
In the approach proposed, the fracture width, radius, and strike are required to be assumed. Incorrect assumptions of fracture width and radius will yield erroneous results, for WA parameters and fracture density, in practice. A priori knowledge about fracture width is important in the prediction of fracture density. However, the interpretation of fracture reservoirs shows that the fracture apertures in reservoir depth are about 1 mm, and have little variation.  where g = V S 2 V p 2 It should be noted that sensitivity F is only related to the elastic properties of the isotropic background.