An Accurate Jacobian Matrix with Exact Zoeppritz for Elastic Moduli of Dry Rock

Abstract

Seismic velocities are related to the solid matrices and the pore fluids. The bulk and shear moduli of dry rock are the primary parameters to characterize solid matrices. Amplitude variation with offset (AVO) or amplitude variation with incidence angle (AVA) is the most used inversion method to discriminate lithology in hydrocarbon reservoirs. The bulk and shear moduli of dry rock, however, cannot be inverted directly using seismic data and the conventional AVO/AVA inversions. The most important step to accurately invert these dry rock parameters is to derive the Jacobian matrix. The combination of exact Zoeppritz and Biot–Gassmann equations makes it possible to directly calculate the partial derivatives of seismic reflectivities (PP-and PS-waves) with respect to dry rock moduli. During this research, we successfully derive the accurate partial derivatives of the exact Zoeppritz equations with respect to bulk and shear moduli of dry rock. The characteristics of these partial derivatives are investigated in the numerical examples. Additionally, we compare the partial derivatives using this proposed algorithm with the classical Shuey and Aki–Richards approximations. The results show that this derived Jacobian matrix is more accurate and versatile. It can be used further in the conventional AVO/AVA inversions to invert bulk and shear moduli of dry rock directly.

1. Introduction

The amplitude variation with offset (AVO) or amplitude variation with incidence angle (AVA) inverts from the measured amplitudes of seismic reflectivities to the rock properties [1,2,3]. Due to the complexity of Zoeppritz equations, most AVO/AVA inversions are based on linear approximations [4,5,6]. The most widely used classical approximations are proposed by Aki–Richards and Shuey [7,8]. These classical approximations, however, are only valid for weak elastic contrasts and relatively small incidence angles [9,10,11]. Actually, there are many target zones that have strong elastic contrast interfaces or long-offset seismic data used. Thus, the accurate partial derivatives of the reflection coefficients with respect to the unknown parameters should be derived for more accurate AVO/AVA inversion [12,13,14].

Rock physics plays a significant role in oil and gas explorations, such as hydrocarbon detection, lithology determination and fluid identification [15,16,17]. Biot–Gassmann equations are widely used in rock physics, which establish the relationship between elastic moduli (bulk and shear moduli) and seismic P-and S-wave velocities [18,19,20,21]. Fluid substitution is a technique used to predict P-, S-wave velocities and density based on the analyses of their porosity and dry rock moduli [22,23,24,25]. The elastic moduli of dry rock could be obtained from the rock physics experiments, well logging or numerical modeling [26,27]. Benson and Wu predicted dry rock bulk modulus, rigidity modulus and reflection coefficients with the application of laboratory rock samples and well logs [28]. Li and Zhang derived the analytical approximations of bulk and shear moduli of dry rock based on the differential effective medium theory [29]. However, the direct inversion of dry rock moduli from seismic data is always a challenging problem.

During recent decades, some researchers proposed to invert the fluid or elastic moduli from seismic data. Goodway et al. provided an AVO approximate equation through incorporating the dry and saturated properties of rock into the P-wave velocity expression, which is named as the lambda-mu-rho (LMR) method [30]. Russell et al. combined this concept with the Biot–Gassmann theory to derive an approximation equation for AVO inversion [31]. Based on Russell’s fluid factor, Yin et al. carried out the Bayesian-based elastic impedance inversion, which is more suitable for seismic data with large incidence angles [32]. Zong et al. further proposed FMR–AVA inversion (Fluid Factor, Mu (Shear modulus), Rho (Density)–Amplitude Variation with Angle) using the approximations [33]. We propose to combine the exact Zoeppritz equations with the Biot–Gassmann equations.

The best way to get an accurate elastic moduli of dry rock from seismic data is to use AVO/AVA inversion with exact Zeoppritz equations. The most difficult part of AVO/AVA inversion based on the exact Zoeppritz equations, however, is to obtain the accurate Jacobian matrix (partial derivatives of the reflection coefficients with respect to dry rock moduli) [34,35]. An accurate algorithm was proposed to compute the gradients of the seismic reflection coefficients with respect to P-, S-wave velocities and density by solving the partial derivative equations of the exact Zoeppritz equations [36,37]. Lu et al. carried out joint PP and PS AVA seismic inversion using exact Zoeppritz equations [38]. Nevertheless, the direct inversion of bulk and shear moduli of dry rock was not taken into consideration in the previous studies. We propose to combine the exact Zoeppritz equations with the Biot–Gassmann equations and derive the accurate Jacobian matrix for inverting dry rock moduli directly without any simplifications. This derived matrix can be used further to invert bulk and shear moduli of dry rock. Due to the usage of exact Zoeppritz equations, it overcomes the weaknesses that the approximate Zoeppritz equations have.

We derive the accurate partial derivatives of reflection coefficients with respect to dry rock parameters (bulk and shear moduli) through combining the exact Zoeppritz equations and the Biot–Gassmann equations. During the numerical tests, the characteristics of these partial derivatives are observed. Additionally, we compare the partial derivatives and their errors using this proposed method with those from the classical Aki–Richards and Shuey approximations.

2. Zoeppritz and Biot–Gassmann Equations

Zoeppritz equations describe the relationship between the seismic reflection coefficients and the elastic rock properties [39,40]. The Zoeppritz equations for only P-wave incidence are written as:

$$AR=B$$

(1)

where

$$R=\left[\begin{array}{l}{R}_{pp}\\ R_{ps}\\ T_{pp}\\ T_{ps}\end{array}\right],B=\left[\begin{array}{l}{b}_1\\ b_2\\ b_3\\ b_4\end{array}\right], \mathrm{and} A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& a_{14}\\ a_{21}& a_{22}& a_{23}& a_{24}\\ a_{31}& a_{32}& a_{33}& a_{34}\\ a_{41}& a_{42}& a_{43}& a_{44}\end{array}\right]$$

(2)

where $R_{pp}$, $R_{ps}$, $T_{pp}$ and $T_{ps}$ are the reflection and transmission coefficients of P-and SV-wave (P–SV conversion), respectively. The elements a and b in the matrix A and B are

$$a_{11}=\sin \alpha , a_{12}=\cos \beta , a_{13}= -\sin {\alpha }^{\prime }, a_{14}=\cos {\beta }^{\prime };$$

$$a_{21}=\cos \alpha , a_{22}= -\sin \beta , a_{23}=\cos {\alpha }^{\prime }, a_{24}=\sin {\beta }^{\prime };$$

$$a_{31}=\cos 2\beta , a_{32}= -\frac{v_{s1}}{v_{p1}}\sin 2\beta , a_{33}= -\frac{{\rho }_{e2}}{{\rho }_{e1}}\frac{v_{p2}}{v_{p1}}\cos 2{\beta }^{\prime }, a_{34}= -\frac{{\rho }_{e2}}{{\rho }_{e1}}\frac{v_{s2}}{v_{p1}}\sin 2{\beta }^{\prime };$$

$$a_{41}=\frac{v_{s1}^2}{v_{p1}^2}\sin 2\alpha , a_{42}=\frac{v_{s1}}{v_{p1}}\cos 2\beta , a_{43}=\frac{{\rho }_{e2}}{{\rho }_{e1}}\frac{v_{s2}^2}{v_{p2}^{}{v}_{p1}}\sin 2{\alpha }^{\prime }, a_{44}= -\frac{{\rho }_{e2}}{{\rho }_{e1}}\frac{v_{s2}}{v_{p1}}\cos 2{\beta }^{\prime };$$

$$b_1= -\sin \alpha , b_2=\cos \alpha , b_3= -\cos 2\beta , b_4=\frac{v_{s1}^2}{v_{p1}^2}\sin 2\alpha .$$

Here, $\alpha$ denotes the incidence and reflection angle of P-wave; $\beta$ is the reflection angle of SV-wave; ${\alpha }^{\prime }$ and ${\beta }^{\prime }$ denote the transmission angles of the P-and SV-waves, respectively; $v_{p1}$, $v_{s1}$, ${\rho }_1$ and $v_{p2}$, $v_{s2}$, ${\rho }_2$ are the P-, S-wave velocities and density in layer 1 and 2, respectively (shown in Figure 1).

The Biot–Gassmann equations describe the relationship between velocity, density and elastic moduli in porous media. It can be expressed as Hilterman, [41]:

$$v_{pi}^2=\left[\left(K_{di}+4/3{\mu }_{di}\right)+\frac{{\left(1-K_{di}/K_{si}\right)}^2}{\left(1-{\varphi }_i-K_{di}/K_{si}\right)/K_{si}+{\varphi }_i/K_{fi}}\right]/{\rho }_{ei}$$

(3)

$$v_{si}^2={\mu }_{di}/{\rho }_{ei}, {\mu }_{ei}={\mu }_{di}$$

(4)

$${\rho }_{ei}=(1-{\varphi }_i){\rho }_{si}+{\varphi }_i{\rho }_f$$

(5)

where $\varphi$ is porosity; $K_f$, $K_s$ and $K_d$ denote the bulk moduli of effective fluid, solid mineral and dry rock, respectively; ${\rho }_e$, ${\rho }_s$ and ${\rho }_f$ are the densities of saturated medium, solid mineral and pore fluid, respectively; ${\mu }_e$ and ${\mu }_d$ denote the shear moduli of saturated medium and dry frame, respectively.

The approximation of bulk modulus of the effective fluid is as demonstrated by Hilterman, [41]:

$$\frac{1}{K_f}=\frac{S_w}{K_w}+\frac{S_o}{K_o}+\frac{S_g}{K_g}$$

(6)

where $K_o$, $K_g$ and $K_w$ are the bulk moduli of oil, gas and water, respectively. $S_o$, $S_g$ and $S_w$ are saturations when pores are filled with oil, gas and water, respectively.

The effective density can be approximated as

$${\rho }_f=S_w{\rho }_w+S_o{\rho }_o+S_g{\rho }_g$$

(7)

where ${\rho }_o$, ${\rho }_g$ and ${\rho }_w$ are the densities of oil, gas and water, respectively.

3. The Partial Derivatives of P- and S-Wave Velocities with Respect to Dry Rock Moduli

To accurately invert the dry rock parameters, the most important step is to derive the Jacobian matrix (partial derivatives of PP and PS reflection coefficients with respect to dry rock moduli). Letting $m_j$ represent the unknown parameters $K_{d1}$, $K_{d2}$, ${\mu }_{d1}$ and ${\mu }_{d2}$, the partial derivatives of the Zoeppritz equations with respect to $m_j$ can be expressed as Liu et al., [11]

$$A\frac{\partial R}{\partial m_j}=-\frac{\partial A}{\partial m_j}R+\frac{\partial B}{\partial m_j}$$

(8)

Equation (8) shows matrix A depends on the incidence angle. The reflectivity R can be calculated using Zoeppritz Equations (1) and (2). The two remaining parameters needed to be calculated are $\frac{\partial A}{\partial m_j}$ and $\frac{\partial B}{\partial m_j}$. The following sections will derive these two partial derivatives.

To derive $\frac{\partial A}{\partial m_j}$ and $\frac{\partial B}{\partial m_j}$, intermediate partial derivatives of P- and S-wave velocities need to be derived first. Using chain rules and Equation (3), the partial derivatives of P-wave velocity with respect to elastic moduli $K_d$ and ${\mu }_d$ of dry rock are:

$$\frac{\partial v_{p1}}{\partial K_{d1}}=\frac{1}{2v_{p1}{\rho }_{e1}}\left[1-\frac{1}{b_{k1}^2}\left(\frac{2}{K_{s1}}\left(1-\frac{K_{d1}}{K_{s1}}\right)b_{k1}-\frac{1}{K_{s1}^2}{\left(1-\frac{K_{d1}}{K_{s1}}\right)}^2\right)\right]$$

(9)

$$\frac{\partial v_{p1}}{\partial {\mu }_{d1}}=\frac{2}{3v_{p1} {\rho }_{e1}}$$

(10)

$$\frac{\partial v_{p2}}{\partial K_{d2}}=\frac{1}{2v_{p2}{\rho }_{e2}}\left[1-\frac{1}{b_{k2}^2}\left(\frac{2}{K_{s2}}\left(1-\frac{K_{d2}}{K_{s2}}\right)b_{k2}^{}-\frac{1}{K_{s2}^2}{\left(1-\frac{K_{d2}}{K_{s2}}\right)}^2\right)\right]$$

(11)

$$\frac{\partial v_{p2}}{\partial {\mu }_{d2}}=\frac{2}{3v_{p2}{\rho }_{e2}}$$

(12)

where $b_{k1}=\left(1-{\varphi }_1-\frac{K_{d1}}{K_{s1}}\right)\frac{1}{K_{s1}}+\frac{{\varphi }_1}{K_{f1}}$ and $b_{k2}=\left(1-{\varphi }_2-\frac{K_{d2}}{K_{s2}}\right)\frac{1}{K_{s2}}+\frac{{\varphi }_2}{K_{f2}}$.

Using Equations (3) and (4), we have:

$$v_{p1}^2{\rho }_{e1}=\left(K_{d1}+\frac{4}{3}{v}_{s1}^2{\rho }_{e1}\right)+\frac{{\left(1-K_{d1}/K_{s1}\right)}^2}{b_{k1}^{}}$$

(13)

Additionally, the partial derivatives of S-wave velocity with respect to elastic modulus $K_d$ and ${\mu }_d$ of dry rock are:

$$\frac{\partial v_{s1}}{\partial K_{d1}}=\frac{3v_{s1}}{8{\mu }_{d1}^{}}\left[2v_{p1}{\rho }_{e1}\frac{\partial v_{p1}}{\partial K_{d1}}+\frac{1}{b_{k1}^2}\left(\frac{2}{K_{s1}}\left(1-\frac{K_{d1}}{K_{s1}}\right)b_{k1}-\frac{1}{K_{s1}^2}{\left(1-\frac{K_{d1}}{K_{s1}}\right)}^2\right)-1\right]$$

(14)

$$\frac{\partial v_{s1}}{\partial {\mu }_{d1}}=\frac{1}{2v_{s1}{\rho }_{e1}},$$

(15)

$$\frac{\partial v_{s2}}{\partial K_{d2}}=\frac{3v_{s2}}{8{\mu }_{d2}^{}}\left[2v_{p2}{\rho }_{e2}\frac{\partial v_{p2}}{\partial K_{d2}}+\frac{1}{b_{k2}^2}\left(\frac{2}{K_{s2}}\left(1-\frac{K_{d2}}{K_{s2}}\right)b_{k2}-\frac{1}{K_{s2}^2}{\left(1-\frac{K_{d2}}{K_{s2}}\right)}^2\right)-1\right]$$

(16)

$$\frac{\partial v_{s2}}{\partial {\mu }_{d2}}=\frac{1}{2v_{s2}{\rho }_{e2}}$$

(17)

4. Calculation of Jacobian Matrix Elements Using Exact Zeoppritz and Biot–Gassmann Equations

4.1. The Calculation of $\frac{\partial A}{\partial m_j}$

According to the Snell’s law $\frac{v_{p1}}{\sin \alpha }=\frac{v_{s1}}{\sin \beta }=\frac{v_{p2}}{\sin {\alpha }^{\prime }}=\frac{v_{s2}}{\sin {\beta }^{\prime }}$, we can get $\sin \beta =\frac{v_{s1}}{v_{p1}}\sin \alpha$, $\sin {\alpha }^{\prime }=\frac{v_{p2}}{v_{p1}}\sin \alpha$, $\sin {\beta }^{\prime }=\frac{v_{s2}}{v_{p1}}\sin \alpha$, $\cos \beta =\sqrt{1-{\left(\frac{v_{s1}}{v_{p1}}\sin \alpha \right)}^2}$, $\cos {\alpha }^{\prime }=\sqrt{1-{\left(\frac{v_{p2}}{v_{p1}}\sin \alpha \right)}^2}$ and $\cos {\beta }^{\prime }=\sqrt{1-{\left(\frac{v_{s2}}{v_{p1}}\sin \alpha \right)}^2}$. $\cos 2\beta$ and $\sin 2\beta$ can be calculated using trigonometric functions $\cos 2\beta =1-2 {\sin }^2\beta$ and $\sin 2\beta =2\sin \beta \cos \beta$. Similarly, we can obtain $\cos 2{\alpha }^{\prime }$, $\sin 2{\alpha }^{\prime }$, $\cos 2{\beta }^{\prime }$ and $\sin 2{\beta }^{\prime }$, respectively. The partial derivatives of matrix A with respect to each parameter will be derived in the following section.

4.1.1. Partial Derivative of Matrix A with Respect to $K_{d1}$

The trigonometric functions of $\alpha$, $\beta$, ${\alpha }^{\prime }$ and ${\beta }^{\prime }$ are related to the bulk modulus $K_{d1}$ of the medium, letting:

$${\xi }_{kd1}=\frac{\partial }{\partial K_{d1}}\left(\frac{v_{s1}^{}}{v_{p1}^{}}\right)=\frac{1}{v_{p1}^2}\left(v_{p1}\frac{\partial v_{s1}}{\partial K_{d1}}-v_{s1}\frac{\partial v_{p1}}{\partial K_{d1}}\right)=\frac{v_{s1}}{v_{p1}}\left({\tau }_{kd1}^{vs1}-{\tau }_{kd1}^{vp1}\right)$$

(18)

where ${\tau }_{kd1}^{vp1}=\frac{1}{v_{p1}}\frac{\partial v_{p1}}{\partial K_{d1}}$ and ${\tau }_{kd1}^{vs1}=\frac{1}{v_{s1}}\frac{\partial v_{s1}}{\partial K_{d1}}$. Then we have:

$$\frac{\partial \sin \beta }{\partial K_{d1}}=\frac{\partial }{\partial K_{d1}}\left(\frac{v_{s1}^{}}{v_{p1}^{}}\right)\sin \alpha =\frac{v_{p1}}{v_{s1}^{}}\frac{\partial }{\partial K_{d1}}\left(\frac{v_{s1}^{}}{v_{p1}^{}}\right)\sin \beta ={\eta }_{kd1}\sin \beta$$

(19)

where ${\eta }_{kd1}=\frac{v_{p1}}{v_{s1}^{}}\frac{\partial }{\partial K_{d1}}\left(\frac{v_{s1}^{}}{v_{p1}^{}}\right)=\frac{v_{p1}}{v_{s1}^{}}{\xi }_{kd1}$. Since $\alpha$ is independent of $K_{d1}$, the partial derivatives of the trigonometric functions of angle $\alpha$ with respect to $K_{d1}$ are equal to zero. The detailed derivations of other partial derivatives are shown in Appendix A.

Thus, the partial derivative of matrix A with respect to $K_{d1}$ is:

$$\frac{\partial A}{\partial K_{d1}}={\tau }_{kd1}^{vp1}\left[\begin{array}{cccc}0& -\frac{{\eta }_{kd1}}{{\tau }_{kd1}^{vp1}}\tan \beta \sin \beta & \sin {\alpha }^{\prime }& \tan {\beta }^{\prime }\sin {\beta }^{\prime }\\ 0& \frac{{\eta }_{kd1}}{{\tau }_{kd1}^{vp1}}\sin \beta & \tan {\alpha }^{\prime }\sin {\alpha }^{\prime }& -\sin {\beta }^{\prime }\\ 2\frac{{\eta }_{kd1}}{{\tau }_{kd1}^{vp1}}(a_{31}-1)& \frac{{\eta }_{kd1}}{{\tau }_{kd1}^{vp1}}{a}_{32}(2-{\tan }^2\beta )& -\left(3a_{33}+\frac{2{\rho }_2}{{\rho }_1}\frac{v_{p2}}{v_{p1}}\right)& a_{34}({\tan }^2{\beta }^{\prime }-2)\\ a_{41}\left(2\frac{{\eta }_{kd1}}{{\tau }_{kd1}^{vp1}}+1\right)& \frac{{\tau }_{kd1}^{vs1}}{{\tau }_{kd1}^{vp1}}{a}_{42}+\frac{2{\eta }_{kd1}}{{\tau }_{kd1}^{vp1}}(a_{42}-v_{s1})& a_{43}\left({\tan }^2{\alpha }^{\prime }-1\right)& -2(a_{44}+\frac{{\rho }_2}{{\rho }_1}{v}_{s2})\end{array}\right]$$

(20)

4.1.2. Partial Derivative of Matrix A with Respect to $K_{d2}$

Since $\alpha$ and $\beta$ are independent of $K_{d2}$, the partial derivatives of the trigonometric functions of angles $\alpha$ and $\beta$ with respect to $K_{d2}$ are zero. The detailed derivations of all partial derivatives also are shown in Appendix A.

Accordingly, the partial derivative of matrix A with respect to $K_{d2}$ is:

$$\frac{\partial A}{\partial K_{d2}}=\left[\begin{array}{cccc}0& 0& -{\tau }_{kd2}^{vp2}\sin {\alpha }^{\prime }& -{\tau }_{kd2}^{vs2}\tan {\beta }^{\prime }\sin {\beta }^{\prime }\\ 0& 0& -{\tau }_{kd2}^{vp2}\tan {\alpha }^{\prime }\sin {\alpha }^{\prime }& {\tau }_{kd2}^{vs2}\sin {\beta }^{\prime }\\ 0& 0& {\tau }_{kd2}^{vp2}{a}_{33}+2{\tau }_{kd2}^{vs2}\left(a_{33}+\frac{{\rho }_2}{{\rho }_1}\frac{v_{p2}^{}}{v_{p1}}\right)& {\tau }_{kd2}^{vs2}{a}_{34}(2-{\tan }^2{\beta }^{\prime })\\ 0& 0& a_{43}\left(2{\tau }_{kd2}^{vs2}-{\tau }_{kd2}^{vp2}{\tan }^2{\alpha }^{\prime }\right)& (3a_{44}+2\frac{{\rho }_2}{{\rho }_1}{v}_{s2}){\tau }_{kd2}^{vs2}\end{array}\right]$$

(21)

where ${\tau }_{kd2}^{vp2}=\frac{1}{v_{p2}}\frac{\partial v_{p2}}{\partial K_{d2}}$ and ${\tau }_{kd2}^{vs2}=\frac{1}{v_{s2}}\frac{\partial v_{s2}}{\partial K_{d2}}$.

4.1.3. Partial Derivative of Matrix A with Respect to ${\mu }_{d1}$

Since only $\alpha$ is independent of ${\mu }_{d1}$, the partial derivatives of the trigonometric functions of angle $\alpha$ with respect to ${\mu }_{d1}$ are zero. Letting:

$${\eta }_{\mu d1}=\frac{v_{p1}}{v_{s1}^{}}\frac{\partial }{\partial {\mu }_{d1}}\left(\frac{v_{s1}^{}}{v_{p1}^{}}\right)={\xi }_{kd1}\frac{v_{p1}}{v_{s1}^{}}$$

(22)

$${\xi }_{\mu d1}=\frac{\partial }{\partial {\mu }_{d1}}\left(\frac{v_{s1}^{}}{v_{p1}^{}}\right)=\frac{1}{v_{p1}^2}\left(v_{p1}\frac{\partial v_{s1}}{\partial {\mu }_{d1}}-v_{s1}\frac{\partial v_{p1}}{\partial {\mu }_{d1}}\right)=\frac{v_{s1}}{v_{p1}^{}}({\tau }_{\mu d1}^{vs1}-{\tau }_{\mu d1}^{vp1})$$

(23)

where ${\tau }_{\mu d1}^{vp1}=\frac{1}{v_{p1}}\frac{\partial v_{p1}}{\partial {\mu }_{d1}}$ and ${\tau }_{\mu d1}^{vs1}=\frac{1}{v_{s1}}\frac{\partial v_{s1}}{\partial {\mu }_{d1}}$. The partial derivatives of the trigonometric functions of angles $\beta$, $\alpha \prime$ and $\beta \prime$ with respect to ${\mu }_{d1}$ are shown in Appendix B.

Therefore, the partial derivative of matrix A with respect to ${\mu }_{d1}$ is:

$$\frac{\partial A}{\partial {\mu }_{d1}}={\tau }_{\mu d1}^{vp1}\left[\begin{array}{cccc}0& -\frac{{\eta }_{\mu d1}}{{\tau }_{\mu d1}^{vp1}}\tan \beta \sin \beta & \sin {\alpha }^{\prime }& \tan {\beta }^{\prime }\sin {\beta }^{\prime }\\ 0& \frac{{\eta }_{\mu d1}}{{\tau }_{\mu d1}^{vp1}}\sin \beta & \tan {\alpha }^{\prime }\sin {\alpha }^{\prime }& -\sin {\beta }^{\prime }\\ 2\frac{{\eta }_{\mu d1}}{{\tau }_{\mu d1}^{vp1}}(a_{31}-1)& \frac{{\eta }_{\mu d1}}{{\tau }_{\mu d1}^{vp1}}{a}_{32}(2-{\tan }^2\beta )& -\left[3a_{33}+\frac{2{\rho }_2}{{\rho }_1}\frac{v_{p2}}{v_{p1}}\right]& a_{34}({\tan }^2{\beta }^{\prime }-2)\\ a_{41}\left(2\frac{{\eta }_{\mu d1}}{{\tau }_{\mu d1}^{vp1}}+1\right)& \frac{{\tau }_{\mu d1}^{vs1}}{{\tau }_{\mu d1}^{vp1}}{a}_{42}+\frac{2{\eta }_{\mu d1}}{{\tau }_{\mu d1}^{vp1}}(a_{42}-v_{s1})& a_{43}\left({\tan }^2{\alpha }^{\prime }-1\right)& -2(a_{44}+\frac{{\rho }_2}{{\rho }_1}{v}_{s2})\end{array}\right]$$

(24)

4.1.4. Partial Derivative of Matrix A with Respect to ${\mu }_{d2}$

Regarding this, $\alpha$ and $\beta$ are independent of ${\mu }_{d2}$, but ${\alpha }^{\prime }$ and ${\beta }^{\prime }$ are dependent on ${\mu }_{d2}$. The detailed derivations also are shown in Appendix B. Letting ${\tau }_{\mu d2}^{vp2}=\frac{1}{v_{p2}}\frac{\partial v_{p2}}{\partial {\mu }_{d2}}$ and ${\tau }_{\mu d2}^{vs2}=\frac{1}{v_{s2}}\frac{\partial v_{s2}}{\partial {\mu }_{d2}}$, the partial derivative of matrix A with respect to ${\mu }_{d2}$ is:

$$\frac{\partial A}{\partial {\mu }_{d2}}=\left[\begin{array}{cccc}0& 0& -{\tau }_{\mu d2}^{vp2}\sin {\alpha }^{\prime }& -{\tau }_{\mu d2}^{vs2}\tan {\beta }^{\prime }\sin {\beta }^{\prime }\\ 0& 0& -{\tau }_{\mu d2}^{vp2}\tan {\alpha }^{\prime }\sin {\alpha }^{\prime }& {\tau }_{\mu d2}^{vs2}\sin {\beta }^{\prime }\\ 0& 0& {\tau }_{\mu d2}^{vp2}{a}_{33}+2{\tau }_{\mu d2}^{vs2}\left(a_{33}+\frac{{\rho }_2}{{\rho }_1}\frac{v_{p2}^{}}{v_{p1}}\right)& {\tau }_{\mu d2}^{vs2}{a}_{34}(2-{\tan }^2{\beta }^{\prime })\\ 0& 0& a_{43}\left(2{\tau }_{\mu d2}^{vs2}-{\tau }_{\mu d2}^{vp2}{\tan }^2{\alpha }^{\prime }\right)& (3a_{44}+2\frac{{\rho }_2}{{\rho }_1}{v}_{s2}){\tau }_{\mu d2}^{vs2}\end{array}\right]$$

(25)

4.2. The Calculation of $\frac{\partial B}{\partial m_j}$

Matrix B contains four elements and the partial derivatives of matrix B with respect to $K_{d1}$, $K_{d2}$, ${\mu }_{d1}$ and ${\mu }_{d2}$ are

$$\frac{\partial B}{\partial K_{d1}}={\left[\begin{array}{cccc}0,& 0,& 4{\eta }_{kd1}{\sin }^2\beta ,& b_4\left(2{\eta }_{kd1}+{\tau }_{kd1}^{vp1}\right)\end{array}\right]}^T$$

(26)

$$\frac{\partial B}{\partial {\mu }_{d1}}={\left[\begin{array}{cccc}0,& 0,& 4{\eta }_{\mu d1}{\sin }^2\beta ,& b_4\left(2{\eta }_{\mu d1}+{\tau }_{\mu d1}^{vp1}\right)\end{array}\right]}^T$$

(27)

$$\frac{\partial B}{\partial K_{d2}}={\left[\begin{array}{cccc}0,& 0,& 0,& 0\end{array}\right]}^T$$

(28)

$$\frac{\partial B}{\partial {\mu }_{d2}}=\frac{1}{v_{s2}}{\left[\begin{array}{cccc}0,& 0,& 0,& 0\end{array}\right]}^T$$

(29)

where

$$\frac{\partial b_4}{\partial K_{d1}}=\frac{\partial }{\partial K_{d1}}\left(\frac{v_{s1}^2{v}_{p1}^{}}{v_{p1}^2}\right)\sin 2\alpha =b_4\left(2{\eta }_{kd1}+{\tau }_{kd1}^{vp1}\right)$$

(30)

5. Numerical Examples

To observe the characteristics of these partial derivatives with exact Zoeppritz equations, a two-layer model filled with water and oil is built. Regarding the first oil–water interface case, the upper and lower sandstone layers are filled with oil and water, respectively. Concerning the second water–oil interface case, the upper and lower layers are filled with water and oil, respectively. The model parameters are: ${\rho }_w=1000$ ${\mathrm{kg}/\mathrm{m}}^3$, ${\rho }_o=710$ ${\mathrm{kg}/\mathrm{m}}^3$, ${\rho }_g=102$ ${\mathrm{kg}/\mathrm{m}}^3$, ${\rho }_{s1}=2250$ ${\mathrm{kg}/\mathrm{m}}^3$, ${\rho }_{s2}=2350$ ${\mathrm{kg}/\mathrm{m}}^3$; $K_w= 2.0967\times {10}^9$ $\mathrm{Pa}$, $K_o= 1.2382\times {10}^9$ $\mathrm{Pa}$,$K_g= 0.0208\times {10}^9$ $\mathrm{Pa}$, $K_{s1}= 45.2\times {10}^9$ $\mathrm{Pa}$, $K_{s2}= 52.2\times {10}^9$ $\mathrm{Pa}$; ${\varphi }_1=0.2$ and ${\varphi }_2=0.15$ Liu et al., [37].

5.1. Oil–Water Interface

Seismic waves penetrate from the oil-saturated sandstone to the water-saturated sandstone. The water, oil and gas saturations in two sandstone layers are $S_{w1}=0$ $S_{o1}=0.6$, $S_{g1}=0.4$, $S_{w2}=1.0$, $S_{o2}=S_{g2}=0$, respectively. The partial derivative curves of reflectivities $R_{pp}$ and $R_{ps}$ with respect to the elastic moduli $K_{d1}$, $K_{d2}$, ${\mu }_{d1}$ and ${\mu }_{d2}$ of dry rock are complex numbers.

Figure 2 shows the real and imaginary values of partial derivatives of reflectivity $R_{pp}$ with respect to the bulk moduli $K_{d1}$ and $K_{d2}$. One can observe that there is a peak (singular point) at the critical angle (${50}^{\circ }$). It is because the partial derivative matrix $\frac{\partial A}{\partial K_{d1}}$ contains $\tan {\alpha }^{\prime }$. When the incidence angle is larger than the critical angle, the results of the partial derivatives become complex numbers rather than real numbers. Due to the existence of this discontinuous point, big errors will be observed around the critical angle when carrying out seismic inversion of bulk moduli.

Figure 3 shows the real and imaginary curves of partial derivatives of reflectivity $R_{PP}$ with respect to the shear moduli ${\mu }_{d1}$ and ${\mu }_{d2}$, respectively. Different from Figure 2, there is a zero point rather than a singular point at the critical angle. The polarities of $R_{PP}$ with respect to ${\mu }_{d1}$ and ${\mu }_{d2}$ are opposite.

Figure 4 shows the partial derivatives of reflectivity $R_{ps}$ with respect to the bulk moduli $K_{d1}$ and $K_{d2}$, respectively. Figure 4a,b show similar characteristics to Figure 2a,b, but their polarities are different. Additionally, a singular point also is found at the critical angle. The absolute value of the partial derivatives of the PP-wave is almost two times larger than that of the PS-wave. One can observe that the slopes of all curves are nearly zero when the incidence angle is smaller than ${40}^{\circ }$, which means the reflectivity $R_{ps}$ is not sensitive to the bulk moduli and not suitable to be used to invert the bulk moduli.

Figure 5 shows the partial derivatives of reflectivity $R_{ps}$ with respect to the shear moduli ${\mu }_{d1}$ and ${\mu }_{d2}$ of the dry rock, respectively. Different from Figure 3, only one zero point is found here. The variation of the “real” part in Figure 5 is smaller than that in Figure 3, when the incidence angle is larger than the critical angle. This phenomenon could be caused by fluid substitution.

5.2. Water–Oil Interface

A water–oil model is built to further check the efficiency of this algorithm. When seismic waves penetrate from the water-saturated layer to the oil-saturated layer, the partial derivatives of reflectivities $R_{pp}$ and $R_{ps}$ with respect to the elastic moduli $K_{d1}$, $K_{d2}$, ${\mu }_{d1}$ and ${\mu }_{d2}$ of dry rock become real numbers rather than complex numbers.

Figure 6 shows the partial derivatives of reflectivity $R_{pp}$ with respect to the elastic moduli $K_{d1}$, $K_{d2}$, ${\mu }_{d1}$ and ${\mu }_{d2}$, respectively. Different from Figure 2 and Figure 3, there is no singular or zero point at the critical angle. Due to the substitution of pore fluids, the curves of partial derivatives of $R_{pp}$ with respect to the elastic moduli are asymmetric. Comparing with other curves, the variation of partial derivative of $R_{pp}$ with respect to ${\mu }_{d2}$ is much smaller (Figure 6a). The enlarged details can be observed in Figure 6b by zooming in the curve.

Figure 7 shows the partial derivatives of the reflectivity $R_{ps}$ with respect to the elastic moduli $K_{d1}$, $K_{d2}$, ${\mu }_{d1}$ and ${\mu }_{d2}$, respectively. Similar to Figure 6, these curves are asymmetric. The partial derivative of $R_{pp}$ with respect ${\mu }_{d1}$ is much larger than others. The enlarged partial derivatives of $R_{pp}$, with respect to $K_{d1}$, $K_{d2}$ and ${\mu }_{d2}$, are shown in Figure 7b.

5.3. Gradient Comparison between the Exact and the Classical Approximations of Zoeppritz Equations

To test the accuracy of this proposed method, we compare the partial derivatives calculated using the exact Zoeppritz equations with those using two classical approximations: Aki–Richards and Shuey [6,7]. Here, the same oil–water model is used. These two Zoeppritz approximations are used to calculate the partial derivatives of the seismic reflectivities with respect to the dry rock properties directly.

Figure 8 compares the partial derivatives of reflectivity $R_{PP}$ with respect to $K_{d1}$ and $K_{d2}$ using different algorithms. $R_{kd1,Zoep}$ and $R_{kd2,Zoep}$ are the results of the exact Zoeppritz equations; $R_{kd1,Aki}$ and $R_{kd2,Aki}$ are the results of the Aki–Richards approximation; $R_{kd1,Shu}$ and $R_{kd2,Shu}$ are the results of the Shuey approximation. Seen in the AVO/AVA inversion, the incidence angles of seismic data are normally smaller than the critical angle. The largest incidence angle is ${50}^{\circ }$ in Figure 8b,d. When the incidence angle is smaller than ${15}^{\circ }$, the values calculated using these three algorithms are very similar. One can observe that differences between the exact solution and two approximations are becoming larger and larger with the increase in the incidence angle. The results prove that the approximations are only suitable for the AVO inversion with small incidence angles.

Figure 9 shows the comparison of the absolute and relative errors of partial derivatives of reflectivity $R_{pp}$ with respect to $K_{d1}$ and $K_{d2}$ using Aki–Richards and Shuey approximations. $Re{R}_{kd1,Aki}$, $Re{R}_{kd2,Aki}$, $E{R}_{kd1,Aki}$ and $E{R}_{kd2,Aki}$ are the relative and absolute errors of the Aki–Richards approximation, respectively. $Re{R}_{kd1,Shu}$, $Re{R}_{kd2,Shu}$, $E{R}_{kd1,Shu}$ and $E{R}_{kd2,Shu}$ are the relative and absolute errors of the Shuey approximation, respectively. Two zero points can be observed when the incidence angles are $0^{\circ }$ and ${90}^{\circ }$. Therefore, there are two discontinuity points on the relative errors which make other curves close to zero. Considering this circumstance, the absolute error curves are more helpful to observe the accuracy of each algorithm. Figure 9b,d show the relative and absolute errors of two classic approximations when the incidence angle is smaller than the critical angle and the errors are still very large.

Figure 10 compares the partial derivatives of $R_{PP}$ with respect to ${\mu }_{d1}$ and ${\mu }_{d2}$ using different algorithms. $R_{\mu 1,Zoep}$ and $R_{\mu 2,Zoep}$ are the results of the exact Zoeppritz equations; $R_{\mu 1,Aki}$ and $R_{\mu 2,Aki}$ are the results of the Aki–Richards approximation; $R_{\mu 1,Shu}$ and $R_{\mu 2,Shu}$ are the results of the Shuey approximation. The enlarged details of partial derivatives of the exact Zeoppritz are shown in Figure 10b. The values of the approximations are still different with the exact Zoeppritz, even when the incidence angle is pretty small.

Figure 11 compares the absolute and relative errors of partial derivatives of reflectivity $R_{pp}$ with respect to ${\mu }_1$ and ${\mu }_2$ using the Aki–Richards and Shuey approximations. $Re{R}_{\mu 1,Aki}$, $Re{R}_{\mu 2,Aki}$, $E{R}_{\mu 1,Aki}$ and $E{R}_{\mu 2,Aki}$ are the relative and absolute errors of the Aki–Richards approximation, respectively. $Re{R}_{\mu 1,Shu}$, $Re{R}_{\mu 2,Shu}$, $E{R}_{\mu 1,Shu}$ and $E{R}_{\mu 2,Shu}$ are the relative and absolute errors of the Shuey approximation, respectively. Figure 11a shows there are three singular points when the incidence angles are $0^{\circ }$, ${50}^{\circ }$ and ${90}^{\circ }$, respectively. Figure 11b,c show the results without these three singular points. The errors of the Shuey and Aki–Richards approximations are still very large when the incidence angle is between $8^{\circ }$ and ${50}^{\circ }$.

6. Discussion

Recently, AVA (AVO) or elastic impedance inversion has been carried out using approximations of the exact Zoeppritz equations. These approximations are only applicable for the formations with weak reflection interfaces and seismic data with small incidence angles (or small offsets). To overcome the weaknesses of these approximations, it is necessary to derive the accurate Jacobian matrix, which consists of the partial derivatives of reflection coefficients with respect to unknown parameters.

Regarding reservoir geophysics, an important step is to establish a theoretical model between the seismic attributes and the physical properties of rocks. Gassmann proposed equations to represent the influences of rock pore and fluid properties on seismic responses and attributes [21]. Although there are many assumptions, the Biot–Gassmann equations are still the most practical and commonly used methods.

Previous research demonstrated these dry rock moduli are mainly obtained indirectly from rock physics experiments, well logging or numerical modeling. Researchers had not built a direct inversion relationship between the theoretical rock physics model and the Zoeppritz equations. Here, we combined the Biot–Gassmann and the Zoeppritz equations to establish the relationship between the seismic wave responses of lithology and dry rock properties. The derived accurate Jacobian matrix, which consists of the partial derivatives of reflection coefficients with respect to bulk and shear moduli of dry rock, can be used in the further AVO/AVA inversion. The characteristics of the partial derivatives were analyzed and indicate that the reflectivity $R_{ps}$ should not be used to invert the bulk moduli. Our future research will incorporate these partial derivatives to the AVA inversion and get the dry rock moduli directly from seismic data.

7. Conclusions

Here, the Biot–Gassmann equations and the exact Zoeppritz equations were combined to calculate the partial derivatives of seismic wave reflection coefficients with respect to the dry rock moduli (bulk and shear moduli). It was an accurate way to get the Jacobian matrix rather than using the Zoeppritz approximations. This accurate matrix was further used in the AVO/AVA inversions to invert bulk and shear moduli of dry rock directly. The numerical tests above indicated that, when a seismic wave refracts from the low impedance medium (rock filled with oil) to the high impedance medium (rock filled with water), the partial derivatives of reflection coefficients, with respect to the elastic moduli of dry rock, had a singular or zero point at the critical angle. When a seismic wave refracted from the high impedance medium (rock filled with water) to the low impedance medium (rock filled with oil), the partial derivatives of reflection coefficients, with respect to the elastic moduli of dry rock, were no longer complex numbers. There was no singular or zero point at the critical angle in this case. Additionally, the reflectivity of the PS-wave was not sensitive to the bulk moduli and not suitable to be used in the inversion of them. The errors of classical Shuey and Aki–Richards approximations were still very large, even when the incidence angles were relatively small. Our proposed algorithm provides an important theoretical basis for direct and accurate inversion of dry rock parameters from seismic data in future AVO/AVA inversions.