TY - JOUR
AU - Greenhalgh,, Mark
AB - SUMMARY To simulate seismic wave propagation in natural porous media, the second-order Biot wave equations in the frequency domain are numerically solved with the (2-D) mixed-staggered grid nine-point stencil finite-difference method using both homogeneous and heterogeneous media approaches. A comparison of numerical results on three synthetic model examples shows that a significant discrepancy may be caused with the homogeneous formulation that ignores the spatial gradients of the medium properties. Such an approach, whilst previously applied, is not suitable to use for reflected/refracted wave investigations. A new and easier method is introduced to obtain the spatial derivative terms in the governing partial differential equation through a rotated coordinate system without using the derivative chain rule. This method is to directly rotate the vectors of the wave equations in the rotated coordinate system to the vectors in the classic coordinate system to constitute the wave equation in the mixed grid system. By using this new method, the exact same results are obtained as those with the classic method using the derivative chain rule. The calculated wavefields of solid particle velocity and fluid flux from the second-order Biot wave equation are used to calculate the pore (fluid) pressure through Biot's first-order constitutive equation. The fluid pressure wavefields are almost entirely devoid of all identified shear waves. Attention is also drawn to the fact that the current finite-difference, frequency-domain modelling methods, designed for the highest frequency and the slowest shear wave velocity in the model, could result in undersampling of the Biot slow P wave. Numerical solutions, Computational seismology, Wave propagation 1 INTRODUCTION The well-known Biot theory (Biot 1956a,b; 1962) has long been accepted as the most popular approach for numerical modelling of seismic wave propagation in porous media. Using numerical algorithms such as the finite-difference method or the finite-element method, Biot's equations can be solved in either the time domain or the frequency domain. Most publications on seismic wave modelling in poroelastic media have concentrated on solutions in the time domain, for example, Carcione & Quiroga-Goode (1995); Dai et al. (1995); Carcione (1996, 1998); Carcione & Helle (1999) and Liu et al. (2009, 2011). The majority of these time-domain formulations make a low-frequency approximation to simplify the convolution operation when dealing with a frequency-dependent dynamic permeability. As pointed out by Dupuy et al. (2011), the low-frequency assumption and approximation can result in a loss of validity, for example, in near-surface, saturated media where fluid/solid interactions prevail. Masson & Pride (2010) used an explicit time-stepping finite-difference scheme to solve Biot's equations across the entire band of frequencies based on the analytical inverse Fourier transform (Plyushchenkov & Turchaninov 2000) of the frequency-dependent dynamic permeability. Their numerical results show that the stability criterion is strongly dependent on the value of fluid viscosity. The benefits of solving the Biot equations in the frequency domain are manifold. For example, the difficult issue of small time integration steps is avoided. Also the separation of the equations into stiff parts and non-stiff parts seems not necessary (Dupuy et al. 2011; Yang & Mao 2017). Moreover, with a frequency-domain formulation, the convolution operation can be circumvented, which makes it convenient to solve the wave equations even when complex frequency dispersion characteristics exist, such as with a partial (patchy) saturation model (White 1975; Johnson 2001) and a double porosity model (Pride & Berryman 2003a,b; Pride et al. 2004). In these models, in addition to dynamic permeability, all the moduli in the equations are frequency dependent and will result in extra convolution operations when solved in the time domain. Liu et al. (2018) suggested a method to upscale the effective Biot theory (Pride et al. 2004) to poroviscoelastic models in which the generalized fractional Zener model is used to represent the attenuation mechanisms. The relaxation functions of the generalized fractional Zener models have a simple form of frequency dependence, and this method also enables the solution of poroviscoelastic wave propagation in the frequency-space domain over the full frequency range. The numerical solution of the effective Biot equations will be investigated in a companion paper. This extended version of Biot theory to a specific form of double porosity model has very similar equations to those of classic Biot theory (Pride et al. 2004; Liu et al. 2018). The solution developed in this paper can be easily upgraded to an effective Biot theory by considering the effective and frequency-dependent elastic moduli and dynamic permeability appropriate to the effective Biot model. Another advantage of solving the wave equations in the frequency domain is that it can be done by parallel calculations since the discrete wave equations for each frequency are independent from those of other frequencies. Notwithstanding such benefits, only a few publications deal with the numerical solution of poroelastodynamics in the frequency domain. This includes Dupuy et al. (2011) who worked with the finite-element method, Liu et al. (2014 for the BISQ model of Dvorkin & Nur 1993) and Yang & Mao (2017) who employed the finite-difference method. As is well known, when solved in the frequency domain, the discretized seismic wave equations reduce to a matrix equation with the left-hand side equal to the product of an impedance matrix and the unknown wavefield vector. The right-hand side comprises simply the source vector. The typical wavefields to be specified for porous media include the solid particle displacement, the fluid flux, the total stress, and the fluid pressure. In the frequency domain, the elastic equations can be solved using either a first-order or second-order system of partial differential equations (PDEs). The advantages of using the second-order frequency-domain approach for elastic (and viscoelastic) cases has been well explained by Stekl & Pratt (1998). They stated that using the second-order differencing scheme could dramatically reduce computational costs. Please note that in the frequency domain, computational efficiency is the essential consideration in solving elastic wave problems for which only two types of wavefield are involved (i.e. stress and velocity). For Biot equations or poroviscoelastic formulations, the issue of computation efficiency is even more acute than in the elastic case. Using the first-order system, the wavefield vector comprises all four wavefield types which may make the size of the system matrix extremely large. It is this massive matrix equation that controls the computational efficiency of the frequency-domain finite-difference method. This is a possible reason why the second-order system is more popular because only two of the wavefields need to be included in the wavefield vector. Furthermore, second-order finite differencing with standard ordering stencil schemes leads to well-structured, banded diagonal matrices (Stekl & Pratt 1998). Theoretically speaking, the first-order and second-order systems should give the same wave modelling results. However, unlike the first-order system (particle velocity-stress finite-difference method), which is free of the spatial derivatives of the physical parameters (e.g. see Dai et al. 1995; Moczo et al. 2014 for the time-domain formulation), by the method of Stekl & Pratt (1998) the second-order systems require that the spatial derivatives of the physical parameters be incorporated for inhomogeneous media. This is called the heterogeneous approach (Dai et al. 1995; Stekl & Pratt 1998), to be contrasted with the homogeneous approach in which the spatial derivatives of the physical parameters are ignored. The methods used by Liu et al. (2014) and Yang & Mao (2017) for the second-order system fall into the homogeneous approach, not the more correct heterogeneous approach. This can be easily verified by the fact that their second-order systems omit the spatial derivatives of the physical parameters. According to Dai et al. (1995), if Biot's second-order equations for heterogeneous porous media are solved by applying the homogeneous approach, this may sometimes result in large errors. Norton (2009) also reported that ignoring the heterogeneity in the medium can lead to significant errors in the transmitted/scattered wavefields. In terms of the actual finite-difference scheme, Liu et al. (2014) and Yang & Mao (2017) use a 25-point weighted-averaging scheme. One advantage is that this method requires only two and a half gridpoints per shortest wavelength for accurate modelling (Shin & Sohn 1998). Obviously, the use of fewer gridpoints per shortest wavelength implies more efficient computation but this may be at the expense of accuracy. One drawback of the 25-point stencil, as stated by Hustedt et al. (2004), is that structures whose size is smaller than one wavelength and irregular interfaces cannot be properly discretized when using just two and a half points per wavelength. To achieve a reasonable tradeoff between computational accuracy and efficiency, our method uses a mixed-staggered grid approach to discretize the second derivative operator on the classical Cartesian coordinate system (Pratt & Worthington 1990) as well as a rotated system (Saenger et al. 2000; Saenger & Bohlen 2004). The mixed-grid and staggered-grid approach has only a nine-point stencil for the spatial derivatives and requires four gridpoints per minimum wavelength for accurate modelling (Hustedt et al. 2004). It was applied for frequency-domain acoustic wave modelling (Jo et al. 1996; Hustedt et al. 2004; Amini & Javaherian 2012). In this paper, we extend the frequency-domain mixed-grid finite-difference approach for acoustic wave modelling to viscoelastic modelling in the frequency-domain for Biot porous media. To avoid the errors from the homogeneous approach, we solve the second-order Biot PDEs with the heterogeneous approach. But we also calculate the waves with the homogeneous approach for purposes of comparison. As stated by Moczo et al. (2014), explicit calculations of the spatial gradients of the medium properties have only limited accuracy. They circumvented the problem to some extent by an averaging process. But we contend that ignoring these spatial gradients is strictly incorrect and will cause a greater loss of accuracy with the algorithm. Here we use the difference between the two approaches to show the discrepancy caused with the homogeneous approach in ignoring the spatial gradients of the medium properties. The second-order Biot wave equations are derived by eliminating the stress and pressure wavefields from the first-order governing equations. This includes the constitutive equations, the conservation of linear momentum equation, and the Darcy flow law. Unlike solving the first-order wave equations in which all the wavefields are solved simultaneously, we first solve the wave equations for the solid particle velocity and the fluid flux. Then, with these results, the stress and pore pressure wavefields can be obtained by applying the constitutive equations. The pore pressure is a conservative field. Pure shear waves cannot affect the pore pressure. Therefore, the calculation of pore pressure provides an extra check on the correctness of the wave solution because all the shear waves identified in the solid particle velocity and/or fluid flux wavefields must disappear in the pore pressure wavefield. In order to attenuate spurious waves reflected from the artificial model (grid) boundaries, the nearly perfectly matched layer (NPML) technique has been implemented (Hu et al. 2007). We will investigate three 2-D example porous models. The first model is an infinite homogeneous model for which the numerical solution can be compared with the known analytic solution to validate our method. The second model is a two-layer structure for which the P and Swavefronts (e.g. direct waves, reflections, or refractions) can be easily identified. In addition, the discrepancy caused by using the homogeneous approach will be investigated. The third model is a three-layer model involving one complex interface (horst-like structure). This model, which entails laterally varying medium properties, is used to assess the errors that can arise when using the homogenous approach. Based on the first model, an example with a high frequency source is employed to illustrate the problem of inadequate sampling of the slow P wave in the frequency-domain solution. 2 THEORY AND METHOD The Biot constitutive equations (Biot 1962; Pride et al. 2004) in the frequency domain (or for harmonic time dependence |${e^{ - i\omega t}}$|⁠) can be written as \begin{eqnarray*} {\left[ \begin{array}{@{}l@{}} {p_c}\\ {p_f} \end{array} \right] = \frac{1}{{i\omega }}\left[ \begin{array}{@{}l@{}} {K^{sat}}\quad \alpha M\\ \alpha M\quad M \end{array} \right]\left[ \begin{array}{@{}l@{}} \nabla \bullet {{\bf v}}\\ \nabla \bullet {{\bf q}} \end{array} \right]} \end{eqnarray*} (1) \begin{eqnarray*} - i\omega {{{\boldsymbol \sigma }}^D} = \mu \left[ {\nabla {{\bf v}} + {{\left( {\nabla {{\bf v}}} \right)}^T} - {\textstyle{2 \over 3}}\nabla \bullet {{\bf vI}}} \right], \end{eqnarray*} (2) where |${{\bf I}}$| is the unit matrix and the bold dot ‘|$\bullet$|’ denotes the vector dot product. We have the solid particle velocity, |${{\bf v}} = - i\omega {{\bf u}}$|⁠, the fluid flux, |${{\bf q}} = - i\omega {{\bf w}}$|⁠. Quantity |${u_i}$| is the displacement of the solid frame component, |${w_i}$| is the displacement of the fluid relative to the solid (measured in volume per unit area) and the deviatoric stress tensor σD is given by \begin{eqnarray*} \sigma _{ij}^D = {\sigma _{ij}} + {p_c}{\delta _{ij}}, \end{eqnarray*} (3) where the total pressure, |${p_c} = {{ - tr\{ {{\sigma _{ij}}} \}} \mathord{/ {\vphantom {{ - tr\{ {{\sigma _{ij}}} \}} 3}} } 3}$| and |${\sigma _{ij}}$| is the total stress in the medium (including the porous solid frame and the fluid filling the pores) with |${p_f}$| being the fluid pressure in the pores. The stress coefficient |$\alpha$|⁠, the coupling modulus |$M$| between the solid frame and the fluid, the other coefficient, |${\lambda _c}$| and the undrained bulk modulus |${K^{\mathrm{ sat}}}$| (Biot & Willis 1957; Brown & Korringa 1975) are given by the following formulae in terms of bulk modulus and shear modulus of the porous frame, |${K^d}$|⁠, |$\mu$| and those of the solid grains and the pore fluid, |${K^g},{K^f}$| as well as the porosity |$\varphi$|⁠: \begin{eqnarray*} \begin{array}{*{20}{c}} {\alpha = 1 - \displaystyle\frac{{{K^d}}}{{{K^g}}};}& \quad {\displaystyle\frac{1}{M} = \displaystyle\frac{\alpha }{{{K^g}}} + \varphi \left( {\displaystyle\frac{1}{{{K^f}}} - \displaystyle\frac{1}{{{K^g}}}} \right);} \end{array} \end{eqnarray*} (4) \begin{eqnarray*} \begin{array}{*{20}{c}} {{\lambda _c} = {K^{\mathrm{ sat}}} - \displaystyle\frac{2}{3}\mu ;}& \quad {{K^{\mathrm{ sat}}} = {K^d} + {\alpha ^2}M.} \end{array} \end{eqnarray*} (5) The conservation of linear momentum equation in the composite medium is given (in the frequency domain) by \begin{eqnarray*} \nabla \bullet {{{\boldsymbol \sigma }}^D} - \nabla {p_c} + {\rm{i}}\omega \left( {\rho {{\bf v}} + {\rho _f}{{\bf q}}} \right) = - {{{\bf F}}^s}. \end{eqnarray*} (6) Here, |${\rho _f}$| is the density of pore fluid; |$\rho = (1 - \varphi ){\rho _s} + \varphi {\rho _f}$|⁠; |${\rho _s}$| is the density of solid grains; |${{{\bf F}}^s}$| is the body force source term. The dynamic Darcy flow law can be written as \begin{eqnarray*} {{\bf q}} = - \frac{{\kappa \left( \omega \right)}}{\eta }\left( {\nabla {p_f} - {\rm{i}}\omega {\rho _f}{{\bf v}} + {{{\bf F}}^f}} \right), \end{eqnarray*} (7) where |$\eta$| is fluid viscosity; |${{{\bf F}}^f}$| is the body force source for the pore fluid. The dynamic permeability |$\kappa ( \omega )$| can be viewed as |${\kappa _0}$| multiplied by its frequency correction factor (Johnson et al. 1978; Pride et al. 2004), \begin{eqnarray*} \kappa \left( \omega \right) = {{{\kappa _0}} \mathord{\left/ {\vphantom {{{\kappa _0}} {\left( {\sqrt {1 - {{i4\omega } \mathord{\left/ {\vphantom {{i4\omega } {\left( {{\omega _r}{n_J}} \right)}}} \right. } {\left( {{\omega _r}{n_J}} \right)}}} - {{i\omega } \mathord{\left/ {\vphantom {{i\omega } {{\omega _r}}}} \right. } {{\omega _r}}}} \right)}}} \right. } {\left( {\sqrt {1 - {{i4\omega } \mathord{\left/ {\vphantom {{i4\omega } {\left( {{\omega _r}{n_J}} \right)}}} \right. } {\left( {{\omega _r}{n_J}} \right)}}} - {{i\omega } \mathord{\left/ {\vphantom {{i\omega } {{\omega _r}}}} \right. } {{\omega _r}}}} \right)}}, \end{eqnarray*} (8) where |${n_J}$| is suggested by Pride et al. (2004) to have a value of 8. The relaxation frequency |${\omega _r}$| is given by \begin{eqnarray*} {\omega _r} = {{\eta \varphi } \mathord{\left/ {\vphantom {{\eta \varphi } {\left( {{\rho _f}{\kappa _0}T} \right)}}} \right. } {\left( {{\rho _f}{\kappa _0}T} \right)}}, \end{eqnarray*} (9) where |$T$|is tortuosity related to porosity and structure of the pore system. Considering the heterogeneity of the material properties and eliminating the stress terms from eqs (1) to (7) leads to the following second-order heterogeneous Biot wave equations: \begin{eqnarray*} \nabla \bullet \left[ {\mu \nabla {{\bf v}} + \mu {{\left( {\nabla {{\bf v}}} \right)}^T} - {\textstyle{2 \over 3}}\mu \nabla \bullet {{\bf vI}}} \right] + \nabla \left[ {{K^{sat}}\nabla \bullet {{\bf v}} + \alpha M\nabla \bullet {{\bf q}}} \right] + {\omega ^2}\left( {\rho {{\bf v}} + {\rho _f}{{\bf q}}} \right) = {{{\bf f}}^s} \end{eqnarray*} (10) \begin{eqnarray*} \nabla \left[ {\alpha M\nabla \bullet {{\bf v}} + M\nabla \bullet {{\bf q}}} \right] + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{{\bf q}} + {\omega ^2}{\rho _f}{{\bf v}} = {{{\bf f}}^f}. \end{eqnarray*} (11) Here |${{{\bf f}}^s}$| and |${{{\bf f}}^f}$| represent |$- i\omega {{{\bf F}}^s}$| and |$- i\omega {{{\bf F}}^f}$|⁠, respectively, which are body force temporal derivative terms in the frequency domain. With the spatial derivatives of the material properties appearing in eqs (10) and (11), we can consider this to be the heterogeneous (medium) formulation (Stekl & Pratt 1998). By ignoring the spatial derivatives on the material properties, we get the second-order Biot wave equations for a homogeneous medium given below by eqs (12) and (13). This is called the homogeneous formulation to distinguish it from the more general heterogeneous formulation. \begin{eqnarray*} \mu {\nabla ^2}{{\bf v}} + \left( {{\textstyle{1 \over 3}}\mu + {K^{sat}}} \right)\nabla \left( {\nabla \bullet {{\bf v}}} \right) + \alpha M\nabla \left( {\nabla \bullet {{\bf q}}} \right) + {\omega ^2}\left( {\rho {{\bf v}} + {\rho _f}{{\bf q}}} \right) = {{{\bf f}}^s} \end{eqnarray*} (12) \begin{eqnarray*} \alpha M\nabla \left( {\nabla \bullet {{\bf v}}} \right) + M\nabla \left( {\nabla \bullet {{\bf q}}} \right) + {\omega ^2}\left( {{\rho _f}{{\bf v}} + \frac{{i\eta }}{{\omega \kappa \left( \omega \right)}}{{\bf q}}} \right) = {{{\bf f}}^f}. \end{eqnarray*} (13) The Biot wave equations for the homogeneous medium are just a special case of the heterogeneous formulation. Our following investigation will be based on the heterogeneous formulation. But we will also show the solution from the homogeneous formulation for comparison purposes. In the 2-D case, eqs (10) and (11) can be written as \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{l}@{}} {{\omega ^2}\rho {v_x} + {\omega ^2}{\rho _f}{q_x} + \sum\limits_j^4 {{A_{1j}}\left( {{v^{\left( j \right)}}} \right)} = f_x^s}&\quad{{\omega ^2}\rho {v_z} + {\omega ^2}{\rho _f}{q_z} + \sum\limits_j^4 {{A_{2j}}\left( {{v^{\left( j \right)}}} \right)} = f_z^s}\\ {{\omega ^2}{\rho _f}{v_x} + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_x} + \sum\limits_j^4 {{A_{3j}}\left( {{v^{\left( j \right)}}} \right)} = f_x^f}&\quad{{\omega ^2}{\rho _f}v{}_z + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_z} + \sum\limits_j^4 {{A_{4j}}\left( {{v^{\left( j \right)}}} \right)} = f_z^f} \end{array}} \right\}, \end{eqnarray*} (14) where |${v^{( j )}}$| represents the elements of the field quantities |${{\bf v}}$| and |${{\bf q}}$|⁠, in the form of the vector |$\left [ {v_x}\,\,{v_z}\,\,{q_x}\,\,{q_z} \right ]^T$| and the operator terms |${A_{ij}}( {{v_j}} )( {i,j = 1,2,3,4} )$| are given in Appendix A. They are the derivative operations on the field quantity |${v^{( j )}}$| with material properties included, for example \begin{eqnarray*} {A_{11}}\left( {{v_x}} \right) = {\partial _x}\left[ {\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3}} \right){\partial _x}{v_x}} \right] + {\partial _z}\left( {\mu {\partial _z}{v_x}} \right). \end{eqnarray*} (15) We will use the (2-D) mixed-grid nine-point stencil and staggered-grid finite-difference method (Stekl & Pratt 1998; Hustedt et al. 2004) to solve the heterogeneous formulation of the second-order Biot equations in the frequency domain. Here, the mixed grid means a combination of two staggered-grid stencils on the classical Cartesian coordinate system and the 45 deg rotated grid (Saenger et al. 2000; Saenger & Bohlen 2004). For the acoustic wave equations, involving just the pressure wavefield, one only needs to deal with a scalar quantity. The combination of wave equations of two coordinate systems can be conducted through direct weighted summation (Jo et al. 1996; Hustedt et al. 2004; Amini & Javaherian 2012). Unlike the acoustic wave equation, the elastic or poroviscoelastic wave equations, i.e. eq. (14), are vector equations. To combine the classic-grid with the rotated-grid approach, the wavefield components must be aligned in the same direction before the weighted summation is carried out. Saenger et al. (2000),Saenger & Bohlen (2004), and Stekl & Pratt (1998) use the derivative chain rule to change only the directions of derivatives but leave the wave equation and the material parameters (stiffness tensor) unchanged. In other words, by this method the spatial derivatives in the classic coordinate systems are calculated through the spatial derivatives in the rotated coordinate system. But we contend that the medium properties are independent of the coordinate system and therefore the wave equations can be developed directly in the rotated coordinate system. Then with linear rotation equations, the wavefield vectors can be transformed back into the classic coordinate frame for weighted summation with the corresponding components. This method does not need to apply the derivate chain rule. Actually, we apply the two different methods in parallel for numerical calculations and obtain exactly the same result, which validates our method to a certain extent. 3 WAVE EQUATIONS WITH A ROTATED COORDINATE SYSTEM: METHOD 1 The rotated coordinate system |${v{^{\prime}}_x}$| and |${v{^{\prime}}_z}$|⁠, at angle θ to the original horizontal axis, can be expressed in the following way: \begin{eqnarray*} \begin{array}{@{}*{1}{c}@{}} {\left\{ \begin{array}{@{}l@{}} {{v{^{\prime}}}_{x{^{\prime}}}} = {v_x}\cos \theta - {v_z}\sin \theta \\ {{v{^{\prime}}}_{z{^{\prime}}}} = {v_x}\sin \theta + {v_z}\cos \theta \end{array} \right.}\qquad \qquad {{\rm{or}}\,\,{{\bf v{^{\prime}}R}} = {{\bf v}}\,\,\,\,{\rm{and}}\,\,\,{{\bf R}} = \left( {\begin{array}{@{}*{2}{c}@{}} {\cos \theta }&{ - \sin \theta }\\ {\sin \theta }&{\cos \theta } \end{array}} \right)}. \end{array} \end{eqnarray*} (16) Here |${{\bf R}}$| is called the rotation matrix. We denote the mesh spacings along the X-axis and Z-axis of the staggered grid by Δx and Δz, respectively. It follows that |$\cos \theta = {{\Delta x} \mathord{/ {\vphantom {{\Delta x} {\sqrt {{{( {\Delta x} )}^2} + {{( {\Delta z} )}^2}} }}} } {\sqrt {{{( {\Delta x} )}^2} + {{( {\Delta z} )}^2}} }}$| and |$\sin \theta = {{\Delta z} \mathord{/ {\vphantom {{\Delta z} {\sqrt {{{( {\Delta x} )}^2} + {{( {\Delta z} )}^2}} }}} } {\sqrt {{{( {\Delta x} )}^2} + {{( {\Delta z} )}^2}} }}$|⁠. If we set |$\Delta x \ne \Delta z$|⁠, then wave equations in a rotated coordinate system rotated by different angles can be obtained (Shin & Sohn 1998). But, here in this paper we set |$\Delta x = \Delta z$|⁠, |$\theta = {45^0}$| (Fig. 1). Figure 1. View largeDownload slide Illustration of the staggered grid nine-point finite-difference stencil used for the mixed classic and rotated coordinate systems. Figure 1. View largeDownload slide Illustration of the staggered grid nine-point finite-difference stencil used for the mixed classic and rotated coordinate systems. The traditional way to obtain the spatial derivatives in the rotated coordinate frame for elastic waves is to use the derivative chain rule with eq. (16) to replace all the derivatives in eq. (14) of the classic coordinate system by the derivatives in the rotated coordinate system (Stekl & Pratt1998). For example, the derivative |${\partial _x}{v_x}$| in eq. (15) is replaced by \begin{eqnarray*} {\partial _{\bar{x}}}{v_x} = \displaystyle\frac{{\partial {v_x}}}{{\partial x{^{\prime}}}}\displaystyle\frac{{\partial x{^{\prime}}}}{{\partial x}} + \displaystyle\frac{{\partial {v_x}}}{{\partial z{^{\prime}}}}\displaystyle\frac{{\partial z{^{\prime}}}}{{\partial x}} = \cos \theta \displaystyle\frac{{\partial {v_x}}}{{\partial x{^{\prime}}}} + \sin \theta \displaystyle\frac{{\partial {v_x}}}{{\partial z{^{\prime}}}}. \end{eqnarray*} (17) Here |${\partial _{\bar{x}}}$| (or |${\partial _{\bar{z}}}$|⁠) implies that the derivative |${\partial _x}$| (or |${\partial _z}$|⁠) in the classic coordinate system is calculated through the derivatives |${\partial _{x{^{\prime}}}}$| and |${\partial _{z{^{\prime}}}}$| which are along the rotated coordinate axes and obtained by the rotation matrix of eq. (16). By this method, the second-order spatial derivatives involving the rotated coordinate system, e.g. |${\partial _{\bar{z}}}( {\mu {\partial _{\bar{z}}}v} )$| are obtained. They are listed in Appendix B. It is worth repeating that the finite-difference scheme in the rotated coordinate system here for elastic waves is different from that given by Jo et al. (1996) for acoustic (scalar) waves in which the derivatives |${\partial _{x{^{\prime}}}}$| and |${\partial _{z{^{\prime}}}}$| are not used to calculate |${\partial _x}$| and |${\partial _z}$|⁠. Replacing all spatial derivatives |${\partial _x}$| and |${\partial _z}$| in eq. (14) with |${\partial _{\bar{x}}}$| and |${\partial _{\bar{z}}}$|⁠, and applying the equations in Appendix B, we obtain \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{\omega ^2}\rho {v_x} + {\omega ^2}{\rho _f}{q_x} + \sum\limits_j^4 {{{\bar{A}}_{1j}}\left( {{v^{\left( j \right)}}} \right)} = f_x^s}&\quad{{\omega ^2}\rho {v_z} + {\omega ^2}{\rho _f}{q_z} + \sum\limits_j^4 {{{\bar{A}}_{2j}}\left( {{v^{\left( j \right)}}} \right)} = f_z^s}\\ {{\omega ^2}{\rho _f}{v_x} + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_x} + \sum\limits_j^4 {{{\bar{A}}_{3j}}\left( {{v^{\left( j \right)}}} \right)} = f_x^f}&\quad{{\omega ^2}{\rho _f}v{}_z + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_z} + \sum\limits_j^4 {{{\bar{A}}_{4j}}\left( {{v^{\left( j \right)}}} \right)} = f_z^f} \end{array}} \right\}. \end{eqnarray*} (18) The operators |${\bar{A}_{ij}}( {{v^{( j )}}} )( {i,j = 1,\,\,2,\,\,3,\,\,4} )$| are listed in Appendix C. 4 WAVE EQUATIONS WITH A ROTATED COORDINATE SYSTEM: METHOD 2 Now, we introduce a new method to obtain the spatial derivatives in the rotated coordinate system without using the derivative chain rule, hence it is an easier formulation. The new method is to directly rotate the vectors of the wave equations in the rotated coordinate system to the vectors in the classic coordinate system. Because the physical properties of the medium are independent on the coordinate system, the second-order Biot equations can be directly written in the rotated coordinate system as \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{\omega ^2}\rho {{v{^{\prime}}}_{x{^{\prime}}}} + {\omega ^2}{\rho _f}{{q{^{\prime}}}_{x{^{\prime}}}} + \sum\limits_j^4 {{{A{^{\prime}}}_{1j}}\left( {{{v{^{\prime}}}^{\left( j \right)}}} \right)} = {f{^{\prime}}}_{x{^{\prime}}}^s}&\quad{{\omega ^2}\rho {{v{^{\prime}}}_{z{^{\prime}}}} + {\omega ^2}{\rho _f}{{q{^{\prime}}}_{z{^{\prime}}}} + \sum\limits_j^4 {{{A{^{\prime}}}_{2j}}\left( {{{v{^{\prime}}}^{\left( j \right)}}} \right)} = {f{^{\prime}}}_{z{^{\prime}}}^s}\\ {{\omega ^2}{\rho _f}{{v{^{\prime}}}_{x{^{\prime}}}} + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{{q{^{\prime}}}_{x{^{\prime}}}} + \sum\limits_j^4 {{{A{^{\prime}}}_{3j}}\left( {{{v{^{\prime}}}^{\left( j \right)}}} \right)} = {f{^{\prime}}}_{x{^{\prime}}}^f}&\quad{{\omega ^2}{\rho _f}{{v{^{\prime}}}_{z{^{\prime}}}} + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{{q{^{\prime}}}_{z{^{\prime}}}} + \sum\limits_j^4 {{{A{^{\prime}}}_{4j}}\left( {{{v{^{\prime}}}^{\left( j \right)}}} \right)} = {f{^{\prime}}}_{z{^{\prime}}}^f} \end{array}} \right\}. \end{eqnarray*} (19) Here, |${v{^{\prime}}^{( j )}}$| represents the element of |$\left [ {v{^{\prime}}\,_x{{\prime}}}\,\,{v{^{\prime}}\,_z{{\prime}}}\,\,{q{^{\prime}}\,_x{{\prime}}}\,\,{q{^{\prime}}\,_z{{\prime}}} \right ]^T$| which is along the rotated coordinate axes. In the equations of the rotated coordinate system, the operators |$A{^{\prime}}_{ij}( {{{v{^{\prime}}}^{( j )}}} )( {i,j = 1,2,3,4} )$| are similar to the operators |${A_{ij}}( {{v_j}} )$| in eqs (14) (but very different from |${\bar{A}_{ij}}$| in eq. 18). The method to get |${A{^{\prime}}_{ij}}( {{{v{^{\prime}}}^{( j )}}} )$| is very simple. The derivative along coordinates |$x,z$| in |${A_{ij}}( {{v_j}} )$| are replaced with the derivatives along the rotated coordinates |$x{^{\prime}},z{^{\prime}}$|⁠; and the wavefields |${v_k},{q_k}$| are replaced with |${v{^{\prime}}_k}$| and |${q{^{\prime}}_k}$|⁠. But the material parameters are kept the same because the rotation of coordinate axes does not change material characteristics. For example, referring to |${A_{11}}( {{v_x}} )$| in eq. (15), we can write \begin{eqnarray*} {A{^{\prime}}_{11}}\left( {{{v{^{\prime}}}_{x{^{\prime}}}}} \right) = {\partial _{x{^{\prime}}}}\left[ {\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3}} \right){\partial _{x{^{\prime}}}}{{v{^{\prime}}}_{x{^{\prime}}}}} \right] + {\partial _{z{^{\prime}}}}\left( {\mu {\partial _{z{^{\prime}}}}{{v{^{\prime}}}_{x{^{\prime}}}}} \right). \end{eqnarray*} (20) It is worth noting that the vectors (e.g. |${f{^{\prime}}}_{x{^{\prime}}}^s$| along |$x{^{\prime}}$|⁠) in eq. (19) in the rotated coordinate system have different directions than the corresponding vectors (e.g. |$f_x^s$| along |$x$|⁠) in eq. (14) for the classic coordinate system. Therefore, we cannot perform a weighted average directly for eqs (14) and (19) to create the mixed grid finite-difference wave equations. To do this, we need to substitute wavefields |${v{^{\prime}}_{k{^{\prime}}}},{q{^{\prime}}_{x{^{\prime}}}}$| and source terms |${f{^{\prime}}}_i^s,\ {f{^{\prime}}}_i^f({k{^{\prime}}} = {x{^{\prime}}},{z{^{\prime}}})$| in eq. (19) with |${v_k}$|⁠, |${q_k}$| and |$f_k^s$|⁠, |$f_k^f$|(⁠|$k = x,z$|⁠) through the rotation matrix (16). For example, items |${A{^{\prime}}_{11}}( {{{v{^{\prime}}}_x}} )$| and |${f{^{\prime}}}_{x{^{\prime}}}^s$| can be rewritten as \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{l}@{}} {{{A}{^{\prime}}_{11}}\left( {{v_x}\cos \theta - {v_z}\sin \theta } \right)}&{ = {\partial _{x{^{\prime}}}}\left[ {\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3}} \right){\partial _{x{^{\prime}}}}\left( {{v_x}\cos \theta - {v_z}\sin \theta } \right)} \right]}\\ {}&{\quad + {\partial _{z{^{\prime}}}}\left( {\mu {\partial _{z{^{\prime}}}}\left( {{v_x}\cos \theta - {v_z}\sin \theta } \right)} \right)} \end{array}} \right\} \end{eqnarray*} (21) \begin{eqnarray*} {f{^{\prime}}}_{x{^{\prime}}}^s = f_x^s\cos \theta - f_z^s\sin \theta . \end{eqnarray*} (22) Then, we multiply |$\cos \theta$| or |$\sin \theta$| with both sides of each of the four eq. (19) and using the equation, |${\cos ^2}\theta + {\sin ^2}\theta = 1$|⁠, we obtain the new rotated second-order Biot equations: \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{\omega ^2}\rho {v_x} + {\omega ^2}{\rho _f}{q_x} + \sum\limits_j^4 {{{R}{^{\prime}}_{1j}}\left( {{v^{\left( j \right)}}} \right)} = f_x^s}&\quad{{\omega ^2}\rho {v_z} + {\omega ^2}{\rho _f}{q_z} + \sum\limits_j^4 {{{R}{^{\prime}}_{2j}}\left( {{v^{\left( j \right)}}} \right)} = f_z^s}\\ {{\omega ^2}{\rho _f}{v_x} + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_x} + \sum\limits_j^4 {{{R}{^{\prime}}_{3j}}\left( {{v^{\left( j \right)}}} \right)} = f_x^f}&\quad{{\omega ^2}{\rho _f}v{}_z + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_z} + \sum\limits_j^4 {{{R}{^{\prime}}_{4j}}\left( {{v^{\left( j \right)}}} \right)} = f_z^f} \end{array}} \right\}. \end{eqnarray*} (23) Here |${R{^{\prime}}_{ij}}( {{v^{( j )}}} )( {i,j = 1,2,3,4} )$| are given in Appendix D. In eq. (23), the wavefields and the source terms are all in the classic coordinate system; the derivative operations are conducted in the rotated coordinate system. For example, we can write \begin{eqnarray*} {R{^{\prime}}_{11}}\left( {{v_x}} \right) = {\cos ^2}\theta \cdot {A{^{\prime}}_{11}}\left( {{v_x}} \right) + \cos \theta \sin \theta \cdot {A{^{\prime}}_{12}}\left( {{v_x}} \right) + \cos \theta \sin \theta \cdot {A{^{\prime}}_{21}}\left( {{v_x}} \right) + {\sin ^2}\theta \cdot {A{^{\prime}}_{22}}\left( {{v_x}} \right) \end{eqnarray*} (24) and \begin{eqnarray*} {A{^{\prime}}_{11}}\left( {{v_x}} \right) = {\partial _{x{^{\prime}}}}\left[ {\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3}} \right){\partial _{x{^{\prime}}}}{v_x}} \right] + {\partial _{z{^{\prime}}}}\left( {\mu {\partial _{z{^{\prime}}}}{v_x}} \right). \end{eqnarray*} (25) It is important to realize that all |${A{^{\prime}}_{ij}}$| terms are very similar to their counterpart |${A_{ij}}$|of eq. (14) in the classic coordinate system, except the spatial derivatives which are carried out in the rotated coordinate system for which the numerical method will be given in the following sections. Although the two wave equations, i.e. eqs (18) and (23) are derived independently, it is not difficult, with simple algebraic operations, to prove that each |${\bar{A}_{ij}}( {{v^{( j )}}} )$| and |${R{^{\prime}}_{ij}}( {{v^{( j )}}} )$| are equal. 5 WAVE EQUATIONS FOR THE MIXED GRID SYSTEM The mixed system of PDEs in heterogeneous porous media take the form: \begin{eqnarray*} \left. \begin{array}{@{}l@{}} {\omega ^2}\rho {v_x} + {\omega ^2}{\rho _f}{q_x} + a\sum\limits_j^4 {{A_{1j}}\left( {{v^{\left( j \right)}}} \right)} + \left( {1 - a} \right)\sum\limits_j^4 {R{^{\prime}}_{1j}\left( {{v^{\left( j \right)}}} \right)} = f_x^s\\ {\omega ^2}\rho {v_z} + {\omega ^2}{\rho _f}{q_z} + a\sum\limits_j^4 {{A_{2j}}\left( {{v^{\left( j \right)}}} \right)} + \left( {1 - a} \right)\sum\limits_j^4 {R{^{\prime}}_{2j}\left( {{v^{\left( j \right)}}} \right)} = f_z^s\\ {\omega ^2}{\rho _f}{v_x} + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_x} + a\sum\limits_j^4 {{A_{3j}}\left( {{v^{\left( j \right)}}} \right)} + \left( {1 - a} \right)\sum\limits_j^4 {R{^{\prime}}_{3j}\left( {{v^{\left( j \right)}}} \right)} = f_x^f\\ {\omega ^2}{\rho _f}v{}_z + \frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}{q_z} + a\sum\limits_j^4 {{A_{3j}}\left( {{v^{\left( j \right)}}} \right)} + \left( {1 - a} \right)\sum\limits_j^4 {R{^{\prime}}_{4j}\left( {{v^{\left( j \right)}}} \right)} = f_z^f \end{array} \right\}, \end{eqnarray*} (26) where a = 0.5461 (Jo et al. 1996; Hustedt et al. 2004). After solving these equations, the solutions for |${v^{( j )}}$| which represent the elements of |${[ {{v_x}\,\,{v_z}\,\,{q_x}\,\,{q_z}} ]^T}$| can be used to obtain the pore pressure |${p_f}$| or the confining pressure |${p_c}$| through eq. (1). The wavefields in the time domain are then obtained through inverse Fourier transformation. 6 NUMERICAL ALGORITHMS The discretization of the wave eq. (26) involves operations on the spatial derivative terms (i.e. |${A_{ij}}$| and |${R}{^{\prime}}_{ij}$|⁠) and the acceleration terms (i.e. |${\omega ^2}\rho v$|⁠, |${\omega ^2}\rho q$| and |$\frac{{i\omega \eta }}{{\kappa ( \omega )}}q$|⁠). To compute the spatial derivatives with a mixed nine-point staggered grid, we follow Hustedt et al. (2004) and Stekl & Pratt (1998) to develop the derivative operation on the classic and rotated coordinate systems. Fig. 1 shows a centre node at |$( {m,n} )$|and its eight neighbouring nodes as well as the intermediate points (see Hustedt et al. 2004 for details). Both coordinate systems are depicted. Now, let us consider the finite-difference operation at the gridpoint (m,n). Following the method of Kelly et al. (1976), and Hustedt et al. (2004), it is not difficult to obtain the classic and rotated nine-point finite-difference stencils. Let |$\rho ( {x,z} )$| represent a material parameter and |$v( {x,z} )$| represent the wavefield. The classic nine-point finite-difference stencils are given by \begin{eqnarray*} \frac{\partial }{{\partial x}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial x}}} \right] \approx \frac{1}{{{{\left( {\Delta x} \right)}^2}}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2},n} \right)\left[ {v\left( {m + 1,n} \right) - v\left( {m,n} \right)} \right]\\ - \rho \left( {m - {1 \mathord{\left/ {\vphantom {1 {2,n}}} \right. } {2,n}}} \right)\left[ {v\left( {m,n} \right) - v\left( {m - 1,n} \right)} \right] \end{array} \right\} \end{eqnarray*} (27) \begin{eqnarray*} \frac{\partial }{{\partial z}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial z}}} \right] \approx \frac{1}{{{{\left( {\Delta z} \right)}^2}}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m,n + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)\left[ {v\left( {m,n + 1} \right) - v\left( {m,n} \right)} \right]\\ - \rho \left( {m,n - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)\left[ {v\left( {m,n} \right) - v\left( {m,n - 1} \right)} \right] \end{array} \right\} \end{eqnarray*} (28) \begin{eqnarray*} \frac{\partial }{{\partial x}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial z}}} \right] \approx \frac{1}{{4\Delta x\Delta z}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m + 1,n} \right)\left[ {v\left( {m + 1,n + 1} \right) - v\left( {m + 1,n - 1} \right)} \right]\\ - \rho \left( {m - 1,n} \right)\left[ {v\left( {m - 1,n + 1} \right) - v\left( {m - 1,n - 1} \right)} \right] \end{array} \right\} \end{eqnarray*} (29) \begin{eqnarray*} \frac{\partial }{{\partial z}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial x}}} \right] \approx \frac{1}{{4\Delta x\Delta z}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m,n + 1} \right)\left[ {v\left( {m + 1,n + 1} \right) - v\left( {m - 1,n + 1} \right)} \right]\\ - \rho \left( {m,n - 1} \right)\left[ {v\left( {m + 1,n - 1} \right) - v\left( {m - 1,n - 1} \right)} \right] \end{array} \right\}. \end{eqnarray*} (30) The rotated nine-point finite-difference stencils are given by \begin{eqnarray*} \frac{\partial }{{\partial x{^{\prime}}}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] \approx \frac{1}{{\Delta {x^2} + \Delta {z^2}}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2},n - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)\left[ {v\left( {m + 1,n - 1} \right) - v\left( {m,n} \right)} \right]\\ - \rho \left( {m - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2},n + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)\left[ {v\left( {m,n} \right) - v\left( {m - 1,n + 1} \right)} \right] \end{array} \right\} \end{eqnarray*} (31) \begin{eqnarray*} \frac{\partial }{{\partial z{^{\prime}}}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right] \approx \frac{1}{{\Delta {x^2} + \Delta {z^2}}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2},n + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)\left[ {v\left( {m + 1,n + 1} \right) - v\left( {m,n} \right)} \right]\\ - \rho \left( {m - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2},n - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}} \right)\left[ {v\left( {m,n} \right) - v\left( {m - 1,n - 1} \right)} \right] \end{array} \right\} \end{eqnarray*} (32) \begin{eqnarray*} \frac{\partial }{{\partial x{^{\prime}}}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right] \approx \frac{1}{{\Delta {x^2} + \Delta {z^2}}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m + \frac{1}{2},n - \frac{1}{2}} \right)\left[ {v\left( {m + 1,n} \right) - v\left( {m,n - 1} \right)} \right]\\ - \rho \left( {m - \frac{1}{2},n + \frac{1}{2}} \right)\left[ {v\left( {m,n + 1} \right) - v\left( {m - 1,n} \right)} \right] \end{array} \right\} \end{eqnarray*} (33) \begin{eqnarray*} \frac{\partial }{{\partial z{^{\prime}}}}\left[ {\rho \left( {x,z} \right)\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] \approx \frac{1}{{\Delta {x^2} + \Delta {z^2}}}\left\{ \begin{array}{@{}l@{}} \rho \left( {m + \frac{1}{2},n + \frac{1}{2}} \right)\left[ {v\left( {m + 1,n} \right) - v\left( {m,n + 1} \right)} \right]\\ - \rho \left( {m - \frac{1}{2},n - \frac{1}{2}} \right)\left[ {v\left( {m,n - 1} \right) - v\left( {m - 1,n} \right)} \right] \end{array} \right\}. \end{eqnarray*} (34) The material parameters at the staggered gridpoints need to be calculated from neighbouring gridpoints with averaging methods. To obtain stable results, arithmetic averaging is used for the density and the attenuation parameter ρ, for example, \begin{eqnarray*} \rho \left( {m + \frac{1}{2},n - \frac{1}{2}} \right) = \frac{1}{2}\left[ {\rho \left( {m,n} \right) + \rho \left( {m + 1,n - 1} \right)} \right], \end{eqnarray*} (35) and harmonic (reciprocal) averaging is required for the modulus parameters u (Moczo et al. 2002), for example, \begin{eqnarray*} {u^{ - 1}}\left( {m + \frac{1}{2},n + \frac{1}{2}} \right) = \frac{1}{4}\left[ {{u^{ - 1}}\left( {m,n} \right) + {u^{ - 1}}\left( {m + 1,n} \right) + {u^{ - 1}}\left( {m + 1,n + 1} \right) + {u^{ - 1}}\left( {m,n + 1} \right)} \right]. \end{eqnarray*} (36) It is worth noting that the finite-difference stencils mentioned above can be easily used to get the corresponding version for the homogeneous approaches (see eqs 12 and 13) by ignoring the derivatives of the material parameters. For example, eq. (27) becomes \begin{eqnarray*} \rho \left( {m,n} \right)\displaystyle\frac{{{\partial ^2}v}}{{\partial {x^2}}} \approx \rho \left( {m,n} \right)\displaystyle\frac{1}{{{{\left( {\Delta x} \right)}^2}}}\left[ {v\left( {m + 1,n} \right) - 2v\left( {m,n} \right) + v\left( {m - 1,n} \right)} \right]. \end{eqnarray*} (37) The operator on the acceleration terms is approximated by \begin{eqnarray*} \left( {{\omega ^2}\rho ,or\displaystyle\frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}} \right)q\left( {m,n} \right) \Rightarrow \left( {{\omega ^2}\rho ,or\displaystyle\frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}} \right)\left[ {c \cdot q\left( {m,n} \right) + d\left( \begin{array}{@{}l@{}} q\left( {m + 1,n} \right) + q\left( {m - 1,n} \right)\\ + q\left( {m,n + 1} \right) + q\left( {m,n - 1} \right) \end{array} \right)} \right]. \end{eqnarray*} (38) Here, |$c = 0.6248$| and |$d = 0.09381$| (Jo et al. 1996; Hustedt et al. 2004). A similar operation is applied to the particle velocity field |$v$|⁠. Since the calculation of pore pressure provides an extra check on the correctness of the wave solution, the numerical algorithm to calculate the pore pressure is given in Appendix E. 7 NUMERICAL EXAMPLES We employ three 2-D example models to illustrate results. The model sizes are all |$300 \cdot \Delta x$| by |$300 \cdot \Delta z$| (namely, 487.5 m by 487.5, where the grid spacings |$\Delta x$| and |$\Delta z$| are set to be 1.625 m). The first model is an infinite fullspace filled with porous material A and the other two layered models comprise porous media A and porous media B (see Table 1 for material properties). Materials A and B have the same properties except for different relaxed bulk modulus and shear modulus values. These differences result in the different relaxed velocities, for material A: Vp0 = 3072 (m s−1), Vs0 = 1606 (m s−1); for material B: Vp0 = 5404 (m s−1), Vs0 = 3407 (m s−1). Table 1. Material properties of the sample rocks and fluids. Grains and fluid Parameter Water Parameter Grains |${K^f}$|(N m−2) 2.3 × 109 |${K^g}$|(N m−2) 3.9 × 1010 |${\rho ^f}$|(kg m−3) 1000 |${\rho ^g}$|(kg m−3) 2650 |${\eta ^f}$|(N m−2) 0.001 – – Sold frame Parameter Material A Material B |${K^d}$|(N m−2) 7.8 × 109 4*7.8 × 109 |$\mu$| (N m−2) 5.94 × 109 4.5*5.94 × 109 |$\varphi$| 0.21 |$\kappa$|(m2) 1.03 × 10−14 |$T$| |$T = 0.5( {1 + {1 \mathord{/ {\vphantom {1 \varphi }} } \varphi }} )$| Grains and fluid Parameter Water Parameter Grains |${K^f}$|(N m−2) 2.3 × 109 |${K^g}$|(N m−2) 3.9 × 1010 |${\rho ^f}$|(kg m−3) 1000 |${\rho ^g}$|(kg m−3) 2650 |${\eta ^f}$|(N m−2) 0.001 – – Sold frame Parameter Material A Material B |${K^d}$|(N m−2) 7.8 × 109 4*7.8 × 109 |$\mu$| (N m−2) 5.94 × 109 4.5*5.94 × 109 |$\varphi$| 0.21 |$\kappa$|(m2) 1.03 × 10−14 |$T$| |$T = 0.5( {1 + {1 \mathord{/ {\vphantom {1 \varphi }} } \varphi }} )$| View Large Table 1. Material properties of the sample rocks and fluids. Grains and fluid Parameter Water Parameter Grains |${K^f}$|(N m−2) 2.3 × 109 |${K^g}$|(N m−2) 3.9 × 1010 |${\rho ^f}$|(kg m−3) 1000 |${\rho ^g}$|(kg m−3) 2650 |${\eta ^f}$|(N m−2) 0.001 – – Sold frame Parameter Material A Material B |${K^d}$|(N m−2) 7.8 × 109 4*7.8 × 109 |$\mu$| (N m−2) 5.94 × 109 4.5*5.94 × 109 |$\varphi$| 0.21 |$\kappa$|(m2) 1.03 × 10−14 |$T$| |$T = 0.5( {1 + {1 \mathord{/ {\vphantom {1 \varphi }} } \varphi }} )$| Grains and fluid Parameter Water Parameter Grains |${K^f}$|(N m−2) 2.3 × 109 |${K^g}$|(N m−2) 3.9 × 1010 |${\rho ^f}$|(kg m−3) 1000 |${\rho ^g}$|(kg m−3) 2650 |${\eta ^f}$|(N m−2) 0.001 – – Sold frame Parameter Material A Material B |${K^d}$|(N m−2) 7.8 × 109 4*7.8 × 109 |$\mu$| (N m−2) 5.94 × 109 4.5*5.94 × 109 |$\varphi$| 0.21 |$\kappa$|(m2) 1.03 × 10−14 |$T$| |$T = 0.5( {1 + {1 \mathord{/ {\vphantom {1 \varphi }} } \varphi }} )$| View Large Fig. 2 shows the source and receiver positions for all three example models as well as the location of the interface of model 2. Note that there is no interface in model 1. Model 3 includes an additional interface with a complex (mound-like) shape. Figure 2. View largeDownload slide The locations of the source, denoted as S, and the receivers denoted as R1 to R6 for all three example models. Model 1 is an infinite homogenous full space filled with material A (or the first layer). Model 2 is a two-layer model with the interface at depth Z = 302.25 m with the second layer comprising material B. Model 3 has an additional (basal) interface with topography that introduces lateral heterogeneity. Figure 2. View largeDownload slide The locations of the source, denoted as S, and the receivers denoted as R1 to R6 for all three example models. Model 1 is an infinite homogenous full space filled with material A (or the first layer). Model 2 is a two-layer model with the interface at depth Z = 302.25 m with the second layer comprising material B. Model 3 has an additional (basal) interface with topography that introduces lateral heterogeneity. A Ricker wavelet with a dominant frequency of 25 Hz is used as the source pulse (vertical-force source). We set the maximum frequency to be 100 Hz. The duration of this wavelet is about 0.1 s. We set the maximum traveltime |${T_{\max }}$| as 0.7 s. According to Shannon's sampling theorem, the frequency interval is set as |${1 \mathord{/ {\vphantom {1 {{T_{\max }}}}} } {{T_{\max }}}}$|⁠. This frequency interval provides 71 frequency sampling points evenly distributed across the spectrum of the Ricker wavelet. Three of the receivers denoted by R1, R2, and R3 are at the source level and distances of 110.5, 221.0 and 331.5 m from the source. The other three receivers denoted by R4, R5 and R6 are aligned in a dipping direction at distances from the source of 156.2, 312.5 and 468.7 m. According to Hustedt et al. (2004), there should be at least four gridpoints per shortest wavelength. Using the shear wave velocity of material A (which has the minimum value) and the maximum frequency of our source wavelet to find the minimum wavelength, the grid spacing (⁠|$\Delta x$|⁠,|$\Delta z$| set to be 1.625 m) corresponds approximately to 10 gridpoints per shortest wavelength. The wavefields are calculated with the heterogeneous formulation derived in this paper and also results are computed using the previously applied homogenous formulation. For economy of space, snapshots are only shown for the heterogeneous formulation. The difference between solutions derived from the two methods will be shown through the waveforms extracted from the six recorders. The differences will convey how ignoring the spatial gradient of the material parameters will adversely affect the wavefield calculation results. All the solid particle velocity waveforms are normalized by the maximum value of the Vz waveform at the first receiver R1. Similarly, all the fluid flux waveforms are normalized with the maximum of the Qz waveform at the first receiver R1. 7.1 Model 1 The first model is a homogeneous full space of material A (see Fig. 2 without the interface that marks the boundary with the second layer). The point force source denoted as S is directed along the Z-axis and located at the point (⁠|$30 \cdot \Delta x,30 \cdot \Delta z$|⁠) or (48.8,48.8). It is very near to the left and the upper edges of the modelling domain, which is a challenge to subdue the grid boundary artificial reflection problem. Because of this, the grid number (or the width) of the NPML zones is set as 53, larger than the value of 40 that was suggested in Amini & Javaherian (2012). Figs 3(a)–(c) show snapshots of the vertical component Vz of the solid particle velocity, the Qz component of fluid flux and the pore pressure Pf, respectively, at a time of 134.0 ms. From these pictures we can easily identify the direct Pwavefront and the direct S wavefront (denoted as dP and dS) and we find that Vz and Qz have similar wavefield patterns. This similarity can be explained with their linear relationship shown in the decoupling of Biot wavefields for homogeneous porous medium (Dutta & Odé 1983; Liu et al. 2018). Fig. 3(c) shows the snapshot of the pore pressure Pf calculated from the solutions of solid particle velocity and fluid flux. The disappearance of direct Swavefronts in the pore pressure field strongly supports our algorithm and verifies the solution. Figure 3. View largeDownload slide View largeDownload slide View largeDownload slide (a) Snapshot of the Vz component of the solid particle velocity at a time of 134.0 ms in Model 1, a homogeneous porous medium. The P and Swavefronts can be easily identified. (b) Snapshot of the Qz component of the fluid flux at a time of 134.0 ms in Model 1, a homogeneous porous medium. (c) Snapshot of the pore pressure Pf at a time of 134.0 ms in Model 1, a homogeneous porous medium. Pf is calculated from the solutions of the solid particle velocity and the fluid flux. Note the disappearance of S waves in the pore pressure field. Figure 3. View largeDownload slide View largeDownload slide View largeDownload slide (a) Snapshot of the Vz component of the solid particle velocity at a time of 134.0 ms in Model 1, a homogeneous porous medium. The P and Swavefronts can be easily identified. (b) Snapshot of the Qz component of the fluid flux at a time of 134.0 ms in Model 1, a homogeneous porous medium. (c) Snapshot of the pore pressure Pf at a time of 134.0 ms in Model 1, a homogeneous porous medium. Pf is calculated from the solutions of the solid particle velocity and the fluid flux. Note the disappearance of S waves in the pore pressure field. However, it should be noted that the slow P waves that should exist as a static mode around the source location (Carcione & Quiroga-Goode 1995; Liu et al. 2009, 2011) do not appear in Fig. 3(c). To the best of our knowledge, there has been no paper published to solve the Biot wave equations in the frequency domain that addresses this problem at low frequencies (at high frequencies the slow P waves change to almost normal propagative Pwaves, e.g. Dupuy et al. 2011; and Yang & Mao 2017). As we know, the existence of slow P waves makes the Biot wave equations ‘stiff’. Thus, the time step needs to be extremely small to maintain stability in the numerical methods for solving the equations (Atkinson et al. 2009). With the research of Carcione & Quiroga-Goode (1995), when the numerical solution is conducted in the time domain using a low-frequency source, the stiff equations can be solved by a second-order time integration scheme based on the implicit Crank-Nicolson method or by a splitting method. With the splitting method, the Biot wave equations are partitioned into two sets of differential equations, one stiff and the other non-stiff. The stiff equations are solved analytically. But in the frequency domain, although the stiff problem seems to be circumvented, the solution could suffer from the problem of inadequate sampling of the slow P wave. For example, with the frequency going up to high sonic and ultrasonic range, the Biot slow P wave could gradually turn from a diffusive wave into a normal propagative wave with phase velocities of up to several hundred metres per second. It is clear that with such low velocities we have to set the grid spacing to satisfy the minimum number of gridpoints per shortest wavelength that depend on the numerical method employed. However, for the slow P wave at low frequency, the wavelength is extremely small. The spatial interval that is normally set using the slowest S wave velocity and the highest frequency does not ensure that the wavefield of the slow P wave would be adequately sampled. This question will be addressed later with an example using a high frequency source. Nevertheless, at low frequencies, the relative movement between the solid frame and the pore fluid is very small. The classic Biot theory shows that seismic waves in a porous solid behave more like viscoelastic waves than elastic waves. Therefore, a numerical method in the frequency domain is still a very important means to solve the Biot wave equations. Actually, it was stated by Carcione & Quiroga-Goode (1995), who worked in the time domain, that all previous studies had not reported the presence of the slow static mode, and had obtained stable solutions despite the fact that some authors used standard explicit schemes with coarse time steps. A similar statement can be made for previous frequency-domain investigations. It is obviously very important to validate our method by first comparing our numerical solution with the corresponding analytical solution. Normally, the analytical solution is obtained by multiplying the Fourier transformed Green function with the source spectrum (in the frequency domain) or convolving the Green's function with the source waveform or relaxation function (in the time domain). We found just two explicit 2-D analytical solutions for seismic waves available: Burridge (1976) for elastic waves in the time domain and Karpfinger et al. (2005) for the Biot waves in the frequency domain. The analytical solutions to our Model 1 based on the two Green's functions give almost the same results for the solid particle velocity waveform at low frequencies as we mentioned above. Waveforms of the normalized solid frame particle velocity of the z-component Vz at the six receivers are shown in Fig. 4(a) (for R1, R2 and R3) and Fig. 4(b) (for R4, R5 and R6). At each receiver, the three waveforms correspond to the numerical solution with the heterogeneous formulation (denoted as Num_Hete), the numerical solution with the homogeneous formulation (denoted as Num_Homo) and the analytical solution (denoted as Analytic based on Karpfinger et al. 2005), respectively. We find that the numerical solutions with the heterogeneous formulation and that with the homogeneous formulation are exactly consistent as they should be when the medium is uniform. Both numerical solutions are in reasonable agreement with the analytical solution. In Fig. 4(a), a slight discrepancy of the Vz waveform at R1 corresponding to the time of direct Pwave may be due to numerical dispersion. We need to investigate this further in future research. However, for all other traces the differences between the analytic and numerical results are minimal and within the expected range for limited order finite difference operators. Figure 4. View largeDownload slide View largeDownload slide (a) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3. The numerical solutions with the heterogeneous formulation denoted as Num_Hete and with homogeneous formulation denoted as Num_Homo are both in reasonable agreement with the analytical solution denoted as Analytic. (b) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R4, R5 and R6. The numerical solutions with the heterogeneous formulation denoted as Num_Hete and with homogeneous formulation denoted as Num_Homo are both in reasonable agreement with the analytical solution denoted as Analytic. Figure 4. View largeDownload slide View largeDownload slide (a) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3. The numerical solutions with the heterogeneous formulation denoted as Num_Hete and with homogeneous formulation denoted as Num_Homo are both in reasonable agreement with the analytical solution denoted as Analytic. (b) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R4, R5 and R6. The numerical solutions with the heterogeneous formulation denoted as Num_Hete and with homogeneous formulation denoted as Num_Homo are both in reasonable agreement with the analytical solution denoted as Analytic. Although the source is very near to the left and upper edges of the model, there is no significant artificial grid boundary reflection observed, implying an effective NPML absorbing boundary condition. 7.2 Model 2 Model 2 is a two-layer model with the interface at Z = 302.25 m, with material A in the upper layer and material B in the second layer (see Table 1 for the material properties and Fig. 2 for the model configuration). The point force source denoted as S is directed along the Zaxis and located at the point (⁠|$30 \cdot \Delta x$|⁠,|$30 \cdot \Delta z$|⁠) or (48.8,48.8). Figs 5(a)–(c) show snapshots of the Vz component of the solid particle velocity, the component of fluid flux and the pore pressure Pf, respectively, at a time of 153.9.0 ms. From these figures, we can easily identify the direct P and direct Swavefronts, the transmitted Pwave and the reflected P wave (denoted as dP, dS, rP and tP, respectively) as well. The direct S wave has a very strong amplitude and interferes with the reflected Pwave. Thus, the reflected P wavefront is significantly blurred in the two figures (i.e. Figs 5a and b). In Fig. 5(c), the calculated pore pressure filters out the shear wave and clearly shows the reflected Pwave together with the direct Pwave and the transmitted Pwave. Figure 5. View largeDownload slide View largeDownload slide View largeDownload slide (a) Snapshot of the Vz component of the solid particle velocity at a time of 153.9 ms in Model 2, a two-layer porous medium. The PP reflected wave (denoted as rP) interferes with the direct Swave (denoted as dS). Symbols dP and tP refer to direct Pwave and transmitted Pwave, respectively. (b) Snapshot of the Qz component of the fluid flux at a time of 153.9 ms in Model 2, a two-layer porous medium. The PP reflected wave (denoted as rP) interferes with the direct Swave (denoted as dS). Symbols dP and tP refer to the direct Pwave and transmitted Pwave, respectively. (c) Snapshot of the pore pressure Pf at 153.9 ms in Model 2, a two-layer porous media. Pf is calculated from the solutions of the solid particle velocity and the fluid flux. Swaves are filtered out in the pore pressure field and the reflected PP wave (denoted as rP) is recovered from its earlier interference with the direct Swave. Figure 5. View largeDownload slide View largeDownload slide View largeDownload slide (a) Snapshot of the Vz component of the solid particle velocity at a time of 153.9 ms in Model 2, a two-layer porous medium. The PP reflected wave (denoted as rP) interferes with the direct Swave (denoted as dS). Symbols dP and tP refer to direct Pwave and transmitted Pwave, respectively. (b) Snapshot of the Qz component of the fluid flux at a time of 153.9 ms in Model 2, a two-layer porous medium. The PP reflected wave (denoted as rP) interferes with the direct Swave (denoted as dS). Symbols dP and tP refer to the direct Pwave and transmitted Pwave, respectively. (c) Snapshot of the pore pressure Pf at 153.9 ms in Model 2, a two-layer porous media. Pf is calculated from the solutions of the solid particle velocity and the fluid flux. Swaves are filtered out in the pore pressure field and the reflected PP wave (denoted as rP) is recovered from its earlier interference with the direct Swave. Figs 6(a) and (b) show the waveforms of the normalized solid frame particle velocity of the z-component Vz at recorders R1, R2 and R3 and recorders R4, R5 and R6, respectively. At each receiver, the waveforms are calculated with the heterogeneous formulation (denoted as Num_Hete) and with the homogeneous formulation (denoted as Num_Homo). Figure 6. View largeDownload slide View largeDownload slide (a) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes discrepancies on the interface-created waves, e.g. the reflected P waves (PP) and reflected Swaves (SS). (b) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R4, R5 and R6. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes discrepancies on the interface-created waves, e.g. the reflected Pwaves (PP), refracted Swaves (FS) and transmitted Swaves (TS). Figure 6. View largeDownload slide View largeDownload slide (a) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes discrepancies on the interface-created waves, e.g. the reflected P waves (PP) and reflected Swaves (SS). (b) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R4, R5 and R6. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes discrepancies on the interface-created waves, e.g. the reflected Pwaves (PP), refracted Swaves (FS) and transmitted Swaves (TS). From Figs 6(a) and (b), we find that waveform solutions with the homogenous formulation have a recognizable discrepancy with the corresponding solutions obtained using the heterogeneous formulation at all interface-created waves, e.g. reflected Pwave (denoted with PP), reflected Swave (denoted with SS), refracted Swave (denoted as FS) and transmitted Swave (denoted TS). The difference between Num_Homo and Num_Hete waveforms was caused by the existence of the interface and the changing material properties. There is a large jump in the gradient at the interface. However, the spatial derivatives with the homogenous formulation are ignored and such velocity and other material property gradients may cause significant discrepancies on seismic events involving the interface. Due to the presence of the interface, the waveform can be viewed as the superposition of the direct waves and the interface-created waves. Comparing with the direct waves, these interface-created waves (e.g. PP waves) normally have relatively weaker signals. For these weak signals, the discrepancy is significant. Moreover, these weaker signals are critically important in seismic exploration. Therefore, the finite-difference method with the homogenous formulation is not suitable for the reflected and refracted wave investigations. 7.3 Model 3 Model 3 comprises three layers filled with Material B in the upper and lowermost layers and Material A in the middle layer. Model 3 is a three-layer model with the first interface at Z = 139.75 m and the second interface at Z = 325 m. This second interface incorporates a rectangular mound structure with its top at Z = 243.75 m and having a width of 162.5 at its centre, giving rise to lateral as well as vertical heterogeneity. By introducing a more complex interface in the form of a rectangular mound, we introduce lateral as well as vertical variations in the medium properties (see Fig. 7). Therefore, model 3 will show the ability of the algorithm to treat a 2-D heterogeneous medium and provide further investigation on the discrepancy caused by the application of the homogenous formulation. Figure 7. View largeDownload slide Illustration of Model 3 having three layers filled with Material B in the upper and bottom layers and Material A in the middle layer. The first interface is at depth Z = 139.75 m and the second interface, at a depth of Z = 3225 m, apart from the middle section that features a mound of width 162.5 m and having its top at a depth Z = 243.75 m. The mound gives rise to lateral heterogeneity of the model. Figure 7. View largeDownload slide Illustration of Model 3 having three layers filled with Material B in the upper and bottom layers and Material A in the middle layer. The first interface is at depth Z = 139.75 m and the second interface, at a depth of Z = 3225 m, apart from the middle section that features a mound of width 162.5 m and having its top at a depth Z = 243.75 m. The mound gives rise to lateral heterogeneity of the model. Figs 8(a)–(c) show the snapshots of the Vz component of the solid particle velocity, the Qz component of fluid flux and the pore pressure Pf, respectively, at a time 153.9.0 ms. From Fig. 8(a) or (b), we can easily identify the direct Swavefront (denoted as dS) at the upper-right corner, the transmitted Swavefront (denoted as tS) in the middle layer, and the transmitted Pwavefront (denoted as tP) in the lowermost layer. The dS and tS waves are totally filtered out from the calculated pore pressure image in Fig. 8(c). The three snapshots (especially Fig. 8c) clearly and accurately show the geometry of the two interfaces. The reason is that the middle layer has much ‘softer’ Biot material than the other layers, which enables larger pore volumetric deformation and hence the larger amplitude of pore pressure wavefields. Figure 8. View largeDownload slide View largeDownload slide View largeDownload slide (a) Snapshot of the Vz component of the solid particle velocity at a time of 153.9 ms in Model 3, a three-layer porous medium. Three wavefronts are easily identified as direct Swaves, transmitted Swave, and transmitted P wave, denoted as dS, tS and tP, respectively. (b) Snapshot of the Qz component of the fluid flux at a time of 153.9 ms in Model 3, a three-layer porous medium. Note the direct S waves (dS), transmitted S wave (tS) and transmitted Pwave (tP). (c) Snapshot of the pore pressure Pf at a time of 153.9 ms in Model 3, a three-layer porous medium. Pf is calculated from the solutions of solid particle velocity and fluid flux. Note that the direct Swave at the upper-right corner in previous two figures has been filtered out in the pore pressure field. Figure 8. View largeDownload slide View largeDownload slide View largeDownload slide (a) Snapshot of the Vz component of the solid particle velocity at a time of 153.9 ms in Model 3, a three-layer porous medium. Three wavefronts are easily identified as direct Swaves, transmitted Swave, and transmitted P wave, denoted as dS, tS and tP, respectively. (b) Snapshot of the Qz component of the fluid flux at a time of 153.9 ms in Model 3, a three-layer porous medium. Note the direct S waves (dS), transmitted S wave (tS) and transmitted Pwave (tP). (c) Snapshot of the pore pressure Pf at a time of 153.9 ms in Model 3, a three-layer porous medium. Pf is calculated from the solutions of solid particle velocity and fluid flux. Note that the direct Swave at the upper-right corner in previous two figures has been filtered out in the pore pressure field. Figs 9(a) and (b) show the waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3 and receivers R4, R5 and R6, respectively. Comparing with the result of the heterogeneous formulation (Num_Hete), we observe (as expected) that the homogeneous formulation (Num_Homo) leads to discernable discrepancies in all the waveforms at the six receivers, especially the significant discrepancies in the waveform of the receiver near interface R4. Figure 9. View largeDownload slide (a) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes discernable discrepancies on all the waveforms at the three receivers. (b) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R4, R5 and R6. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes significant discrepancies to the waveforms at the receivers near interface R4. Figure 9. View largeDownload slide (a) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R1, R2 and R3. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes discernable discrepancies on all the waveforms at the three receivers. (b) Waveforms of the normalized solid frame particle velocity of the z-component Vz at receivers R4, R5 and R6. Comparing with the result of the heterogeneous formulation (Num_Hete), the homogeneous formulation (Num_Homo) causes significant discrepancies to the waveforms at the receivers near interface R4. 7.4 An example with a high frequency source As mentioned above, with increasing frequency the slow P wave gradually changes from a diffusive disturbance towards a propagative wave. According to the dispersion formulae for the decoupled Biot wave equations (Dutta & Ode 1983; Liu et al. 2018), material A has a phase velocity for the slow P wave of 480 m s−1, and an attenuation factor of 1/Q = 0.55 at a frequency fc = 1.0E6 Hz. This attenuation factor indicates that the slow P wave is still not a completely normal propagative wave. We set the grid spacing (⁠|$\Delta x$|⁠,|$\Delta z$|⁠) to be 0.08 mm corresponding approximately to six gridpoints per slow P wavelength. Fig. 10 shows the snapshot of the Qz component of the fluid flus at a time of 3.6 μs. The slow P wave (denoted as sP) and the direct Swave (denoted as dS) can be easily identified. However, this method is difficult to be used in Model 1 with a source centre frequency of 25 Hz at which the phase velocity of the slow P wave is just 4.91 m s−1, and the attenuation factor is 1/Q = 2.0E + 4. The extremely high attenuation factor indicates that the slow P wave rapidly diffuses out. In spite of this, with the criteria for a propagating wave, i.e. at least four gridpoints per wavelength, the grid spacing should be less than 0.05 m. On the other hand, we still need to adjust the frequency interval to ensure enough frequency sampling points (e.g. 71 points) over the spectrum of the Ricker wavelet. But in Models 1 to 3 we set the grid spacing as 3.25 m that satisfies the normal S or Pwave velocities. Apparently, it is hard to find a grid spacing suitable for normal body waves and slow P wave as well. This problem needs further investigation. Figure 10. View largeDownload slide Snapshot of the Qz component of the fluid flux at a time of 3.60 μs in Model 1, but the source centre frequency is 1.0E6 Hz. The slow P wave (denoted as sP) and the direct Swave (denoted as dS) can be easily identified. Figure 10. View largeDownload slide Snapshot of the Qz component of the fluid flux at a time of 3.60 μs in Model 1, but the source centre frequency is 1.0E6 Hz. The slow P wave (denoted as sP) and the direct Swave (denoted as dS) can be easily identified. 8 CONCLUSIONS The heterogeneous formulation of the second-order Biot equations in the frequency domain is developed for a mixed system (classical Cartesian coordinate system and a rotated coordinate system) to provide a finite-difference solution of the PDEs. Two methods, with and without using the derivative chain rule, are used to derive the wave equations through the rotated coordinate system and yield exactly the same results. The (2-D) staggered mixed-grid nine-point stencil finite difference method is used to solve the second-order wave equations with the heterogeneous formulation (and a homogeneous formulation as well). Three 2-D example porous media synthetic models are employed to illustrate our method. The first model is an infinite uniform full space for which the computed wavefields compare favourably with the analytical solution, thus validating our numerical solution. The other two layered models show that significant discrepancies may be caused with the homogeneous formulation that ignores the spatial gradients of the medium properties. We find that such an approach produces spurious results near interfaces and therefore is unsuitable in dealing with reflected and refracted waves. From the solutions of solid particle velocity and fluid flux, the snapshots of the pore pressure are calculated for all the three sample models. The disappearance of all the identified Swavefronts in the pore pressure field strongly supports the validity of our algorithm. However, the slow P wave that should exist as a static disturbance mode around the source location does not appear in our pore pressure snapshot. This fact raises our concern that current finite-difference methods in the frequency domain do not adequately sample the slow P waves and thus cannot fully recover the slow P waves. An example with a high frequency source is employed to provide further insight on this issue. ACKNOWLEDGEMENTS We are grateful to the College of Petroleum Engineering & Geosciences, King Fahd University of Petroleum and Minerals, Kingdom of Saudi Arabia, for supporting this research. We thank Dr Alexey Stovas and Dr Peter Moczo for constructive and insightful comments that have helped to improve the manuscript. REFERENCES Amini N. , Javaherian A. , 2012 . A MATLAB-based frequency-domain finite-difference package for solving 2D visco-acoustic wave equation , Waves Random Complex Media , 21 ( 1 ), 161 – 183 .. https://doi.org/10.1080/17455030.2010.537708 Google Scholar Crossref Search ADS Atkinson K. , Han W. , Stewart D. , 2009 . Numerical Solution of Ordinary Differential Equations , pp. 127 – 131 ., A John Wiley & Sons, Inc . Biot M.A. , 1956a . Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range , J. acoust. Soc. Am. , 28 ( 2 ), 168 – 178 .. https://doi.org/10.1121/1.1908239 Google Scholar Crossref Search ADS Biot M.A. , 1956b . Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher-frequency range , J. acoust. Soc. Am. 28 ( 2 ), 178 – 191 . Biot M.A. , 1962 . Mechanics of deformation and acoustic propagation in porous media , J. Appl. Phys. , 33 ( 4 ), 1482 – 1498 .. https://doi.org/10.1063/1.1728759 Google Scholar Crossref Search ADS Biot M.A. , Willis D.G. , 1957 . The elastic coefficients of the theory of consolidation , J. Appl. Mech. , 24 , 594 – 601 . Brown R.J.S. , Korringa J. , 1975 . On the dependence of the elastic properties of a porous rock on compressibility of the pore fluid , Geophysics , 40 ( 4 ), 608 – 616 .. https://doi.org/10.1190/1.1440551 Google Scholar Crossref Search ADS Burridge R. , 1976 . Some Mathematical topics in seismology , Courant Institute of Mathematical Sciences , New York University . Carcione J.M. , 1996 . Wave propagation in anisotropic, saturated porous media: plane-wave theory and numerical simulation , J. acoust. Soc. Am. , 99 ( 5 ), 2655 – 2666 .. https://doi.org/10.1121/1.414809 Google Scholar Crossref Search ADS Carcione J.M. , 1998 . Viscoelastic effective rheologies for modelling wave propagation in porous media , Geophys. Prospect. , 46 , 249 – 270 .. https://doi.org/10.1046/j.1365-2478.1998.00087.x Google Scholar Crossref Search ADS Carcione J.M. , Helle H.B. , 1999 . Numerical solution of the poroviscoelastic wave equation on a staggered mesh , J. Comput. Phys. , 154 ( 2 ), 520 – 527 .. https://doi.org/10.1006/jcph.1999.6321 Google Scholar Crossref Search ADS Carcione J.M. , Quiroga-Goode G. , 1995 . Some aspects of the physics and numerical modelling of Biot compressional waves , J. Comput. Acoust. , 3 ( 4 ), 261 – 280 .. https://doi.org/10.1142/S0218396X95000136 Google Scholar Crossref Search ADS Dai N. , Vafidis A. , Kanasewich E.R. , 1995 . Wave propagation in heterogeneous, porous media: a velocity-stress, finite-difference method , Geophysics , 60 ( 2 ), 327 – 340 .. https://doi.org/10.1190/1.1443769 Google Scholar Crossref Search ADS Dovorkin J. , Nur A. , 1993 . Dynamic poroelasticity: a unified model with the squirt and the Biot mechanisms , Geophysics , 58 ( 4 ), 524 – 533 .. https://doi.org/10.1190/1.1443435 Google Scholar Crossref Search ADS Dupuy B. , Barros L.D. , Garambois S. , Virieux J. , 2011 . Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous Galerkin method , Geophysics , 76 ( 4 ), N13 – N28 .. https://doi.org/10.1190/1.3581361 Google Scholar Crossref Search ADS Dutta N.C. , Od´e H. , 1983 . Seismic reflections from a gas–water contact , Geophysics , 48 ( 2 ), 148 – 162 .. https://doi.org/10.1190/1.1441454 Google Scholar Crossref Search ADS Hustedt B. , Operto S. , Virieux J. , 2004 . Mixed-grid and staggered-grid finite-difference methods for frequency-domain acoustic wave modelling , Geophys. J. Int. , 157 , 1269 – 1296 .. https://doi.org/10.1111/j.1365-246X.2004.02289.x Google Scholar Crossref Search ADS Hu W.Y. , Abubakar A. , Habashy T.M. , 2007 . Application of the nearly perfectly matched layer in acoustic wave modeling , Geophysics , 72 ( 5 ), SM169 – SM175 .. https://doi.org/10.1190/1.2738553 Google Scholar Crossref Search ADS Jo C.H. , Shin C.S. , Suh J.H. , 1996 . An optimal 9-point, finite-difference, frequency-space, 2-D scalar wave extrapolator , Geophysics , 61 ( 2 ), 529 – 537 .. https://doi.org/10.1190/1.1443979 Google Scholar Crossref Search ADS Johnson D.L. , 2001 . Theory of frequency dependent acoustics in patchy-saturated porous media , J. acoust. Soc. Am. , 110 , 682 – 694 .. https://doi.org/10.1121/1.1381021 Google Scholar Crossref Search ADS Johnson D.L. , Koplik J. , Dashen R. , 1978 . Theory of dynamic permeability and tortuosity in fluid-saturated porous media , J. Fluid. Mech. , 176 , 379 – 402 .. https://doi.org/10.1017/S0022112087000727 Google Scholar Crossref Search ADS Karpfinger F. , , Muller T.M. , Gurevich B. , 2005 . Radiation patterns of seismic waves in poroelastic media , in SEG/Houston 2005 Annual Meeting 1791, SM P1.4 , SEG , Tulsa , OnePetro , pp. 1791 – 1795 . Kelly K.R. , Ward R.W. , Treitel S. , Alford R.M. , 1976 . Synthetic seismograms: a finite-difference approach , Geophysics , 41 , 2 – 27 . Google Scholar Crossref Search ADS Liu C. , Yang Q.J. , Guo Z.Q. , Liu Y. , Lan H.T. , Geng M.X. , Wang D. , 2014 , Simulation of wave propagation in two-phase porous media using a frequency-space domain method , Chin. J. Geophys. (in Chinese) , 57 ( 9 ), 2885 – 2899 . Liu X. , Greenhalgh S. , Wang Y.H. , 2011 . 2.5D poroviscoseismic wave modeling in double porosity media , Geophys. J. Int. , 186 , 1285 – 1294 .. https://doi.org/10.1111/j.1365-246X.2011.05106.x Google Scholar Crossref Search ADS L , iu X. , Greenhalgh S. , Zhou B. , 2009 . Transient solution for poroviscoacoustic wave propagation in double porosity media and its limitations , Geophys. J. Int. , 178 , 375 – 393 .. https://doi.org/10.1111/j.1365-246X.2009.04144.x Google Scholar Crossref Search ADS Liu X. , Greenhalgh S. , Zhou B. , Greenhalgh M. , 2018 . Effective Biot theory and its generalization to poroviscoelastic methods , Geophys. J. Int. , 212 , 1255 – 1273 .. https://doi.org/10.1093/gji/ggx460 Google Scholar Crossref Search ADS Masson Y. , Pride S. , 2010 . Finite-difference modeling of Biot's poroelastic equations across all frequencies , Geophysics , 75 ( 2 ), N33 – N41 .. https://doi.org/10.1190/1.3332589 Google Scholar Crossref Search ADS Moczo P. , Kristek J. , Galis M. , 2014 . The Finite-Difference Modelling of Earthquake Motions, Waves and Ruptures , Cambridge Univ. Press . Moczo P. , Kristek J. , Vavrycuk V. , Archuleta R. , Halada L. , 2002 . Heterogeneous staggered-grid finite-difference modeling of seismic motion with volume harmonic and arithmetic averaging of elastic moduli and densities , Bull. seism. Soc. Am. , 92 , 3042 – 3066 .. https://doi.org/10.1785/0120010167 Google Scholar Crossref Search ADS Norton G.V. , 2009 . Comparison of homogeneous and heterogeneous modeling of transient scattering from dispersive media directly in the time domain , Math. Comput. Simul , 80 , 682 – 692 .. https://doi.org/10.1016/j.matcom.2009.08.018 Google Scholar Crossref Search ADS Plyushchenkov B. D. , Turchaninov V. , 2000 . Acoustic logging modeling by refined Biot's equations , Int. J. Mod. Phys. C , 12 , 305 – 396 . Pratt R.G. , Worthington M.H. , 1990 . Inverse theory applied to multisource cross-hole tomography, part 1: Acoustic wave-equation method , Geophys. Prospect. , 38 , 287 – 310 .. https://doi.org/10.1111/j.1365-2478.1990.tb01846.x Google Scholar Crossref Search ADS Pride S.R. , Berryman J.G. , 2003a . Linear dynamics of double-porosity dual-permeability materials. I. Governing equations and acoustic attenuation , Phys. Rev. , E68 ( 036603 ), 1 – 10 . Pride S.R. , Berryman J.G. , 2003b . Linear dynamics of double-porosity dual-permeability materials. II. Fluid transport equations , Phys. Rev. , E 68 ( 036604 ), 1 – 10 . Pride S.R. , Berryman J.G. , Harris J.M. , 2004 . Seismic attenuation due to wave-induced flow , J. geophys. Res. , 109 ( B01201 ), 1 – 19 . Saenger E.H. , Bohlen T. , 2004 . Finite-difference modeling of viscoelastic and anisotropic wave propagation using the rotated staggered grid , Geophysics , 69 ( 2 ), 584 – 591 .. https://doi.org/10.1190/1.1707078 Google Scholar Crossref Search ADS Saenger E.H. , Gold N. , Shapiro A. , 2000 . Modeling the propagation of elastic waves using a modified finite-difference grid , Wave Motion , 31 , 77 – 92 .. https://doi.org/10.1016/S0165-2125(99)00023-2 Google Scholar Crossref Search ADS Shin C. , Sohn H. , 1998 . A frequency-space 2-d scalar wave extrapolator using extended 25-point finite-difference operator , Geophysics , 63 , 289 – 296 .. https://doi.org/10.1190/1.1444323 Google Scholar Crossref Search ADS Stekl I , Pratt R.G. , 1998 . Accurate viscoelastic modelling by frequency-domain finite-differences using rotated operators , Geophysics , 63 ( 5 ), 1779 – 1794 .. https://doi.org/10.1190/1.1444472 Google Scholar Crossref Search ADS White J.E. , 1975 . Computed seismic speeds and attenuation in rocks with partial gas saturation , Geophysics , 40 ( 2 ), 224 – 232 .. https://doi.org/10.1190/1.1440520 Google Scholar Crossref Search ADS Yang Q.J. , Mao W.J. , 2017 . Simulation of seismic wave propagation in 2-D poroelastic media using weighted-averaging finite-difference stencils in the frequency-space domain , Geophys. J. Int. , 208 , 148 – 161 .. https://doi.org/10.1093/gji/ggw380 Google Scholar Crossref Search ADS APPENDIX A: ELEMENTS OF THE DERIVATIVE OPERATOR MATRIX The operator terms |${A_{ij}}( {{v_j}} )( {i,j = 1,2,3,4} )$| in eq. (14) are derived from eqs (10) and (11) which, in the 2-D case, can be written as the following four equations, sequentially: \begin{eqnarray*} \left\{ \begin{array}{@{}l@{}} \left[ \begin{array}{@{}l@{}} {\partial _x}\left( {\mu {\partial _x}{v_x}} \right) + {\partial _z}\left( {\mu {\partial _z}{v_x}} \right)\\ {\partial _x}\left( {\mu {\partial _x}{v_z}} \right) + {\partial _z}\left( {\mu {\partial _z}{v_z}} \right) \end{array} \right] + \left[ \begin{array}{@{}l@{}} {\partial _x}\left( {\mu {\partial _x}{v_x}} \right) + {\partial _z}\left( {\mu {\partial _x}{v_z}} \right)\\ {\partial _z}\left( {\mu {\partial _z}{v_x}} \right) + {\partial _z}\left( {\mu {\partial _z}{v_z}} \right) \end{array} \right] - {\textstyle{2 \over 3}}\left[ \begin{array}{@{}l@{}} {\partial _x}\left( {\mu \left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right)} \right)\\ {\partial _z}\left( {\mu \left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right)} \right) \end{array} \right] + \\ \\ \left[ \begin{array}{@{}l@{}} {\partial _x}\left( {{K^{sat}}\left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right)} \right)\\ {\partial _z}\left( {{K^{sat}}\left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right)} \right) \end{array} \right] + \left[ \begin{array}{@{}l@{}} {\partial _x}\left( {\alpha M\left( {{\partial _x}{q_x} + {\partial _z}{q_z}} \right)} \right)\\ {\partial _z}\left( {\alpha M\left( {{\partial _x}{q_x} + {\partial _z}{q_z}} \right)} \right) \end{array} \right] + {\omega ^2}\rho \left[ \begin{array}{@{}l@{}} {v_x}\\ {v_z} \end{array} \right] + {\omega ^2}{\rho _f}\left[ \begin{array}{@{}l@{}} {q_x}\\ {q_z} \end{array} \right] = \left[ \begin{array}{@{}l@{}} f_x^s\\ f_z^s \end{array} \right] \end{array} \right. \end{eqnarray*} (A1) \begin{eqnarray*} \left[ \begin{array}{@{}l@{}} {\partial _x}\left( {\alpha M\left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right)} \right)\\ {\partial _z}\left( {\alpha M\left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right)} \right) \end{array} \right] + \left[ \begin{array}{@{}l@{}} {\partial _x}\left( {M\left( {{\partial _x}{q_x} + {\partial _z}{q_z}} \right)} \right)\\ {\partial _z}\left( {M\left( {{\partial _x}{q_x} + {\partial _z}{q_z}} \right)} \right) \end{array} \right] + {\omega ^2}{\rho _f}\left[ \begin{array}{@{}l@{}} {v_x}\\ {v_z} \end{array} \right] + \displaystyle\frac{{i\omega \eta }}{{\kappa \left( \omega \right)}}\left[ \begin{array}{@{}l@{}} {q_x}\\ {q_z} \end{array} \right] = \left[ \begin{array}{@{}l@{}} f_x^f\\ f_z^f \end{array} \right]. \end{eqnarray*} (A2) For the derivative operator |${A_{ij}}( {{v_j}} )( {i,j = 1,2,3,4} )$|⁠, the first subscript implies the |${i^{th}}$| equation of the equations (A1) and (A2), and the second refers to the |${j^{th}}$| field quantity in the vector |${[ {{v_x}\,\,{v_z}\,\,{q_x}\,\,{q_z}} ]^T}$|⁠. For example, the operator |${A_{11}}( {{v_x}} )$| consists of all the derivative terms on the |${v_x}$| in the first equation. \begin{eqnarray*} \left\{ {\begin{array}{@{}*{1}{l}@{}} {{A_{11}}\left( {{v_x}} \right) = {\partial _x}\left[ {\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3}} \right){\partial _x}{v_x}} \right] + {\partial _z}\left( {\mu {\partial _z}{v_x}} \right)}\\ {{A_{12}}\left( {{v_z}} \right) = {\partial _x}\left[ {\left( {{K^{sat}} - {\textstyle{2 \over 3}}\mu } \right){\partial _z}{v_z}} \right] + \left[ {{\partial _z}\left( {\mu {\partial _x}{v_z}} \right)} \right]}\\ {{A_{13}}\left( {{q_x}} \right) = {\partial _x}\left[ {\alpha M{\partial _x}{q_x}} \right]}\\ {{A_{14}}\left( {{q_z}} \right) = {\partial _x}\left[ {\alpha M{\partial _z}{q_z}} \right]} \end{array}} \right. \end{eqnarray*} (A3) \begin{eqnarray*} \left\{ {\begin{array}{@{}*{1}{l}@{}} {{A_{21}}\left( {{v_x}} \right) = {\partial _x}\left( {\mu {\partial _z}{v_x}} \right) + {\partial _z}\left[ {\left( {{K^{sat}} - {{2\mu } \mathord{\left/ {\vphantom {{2\mu } 3}} \right. } 3}} \right){\partial _x}{v_x}} \right]}\\ {{A_{22}}\left( {{v_z}} \right) = {\partial _x}\left( {\mu {\partial _x}{v_z}} \right) + {\partial _z}\left[ {\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3}} \right){\partial _z}{v_z}} \right]}\\ {{A_{23}}\left( {{q_x}} \right) = {\partial _z}\left( {\alpha M{\partial _x}{q_x}} \right)}\\ {{A_{24}}\left( {{q_z}} \right) = {\partial _z}\left( {\alpha M{\partial _z}{q_z}} \right)} \end{array}} \right. \end{eqnarray*} (A4) \begin{eqnarray*} \left\{ \begin{array}{@{}l@{}} {A_{31}}\left( {{v_x}} \right) = {\partial _x}\left[ {\alpha M{\partial _x}{v_x}} \right]\\ {A_{32}}\left( {{v_z}} \right) = {\partial _x}\left[ {\alpha M{\partial _z}{v_z}} \right]\\ {A_{33}}\left( {{q_x}} \right) = {\partial _x}\left[ {M{\partial _x}{q_x}} \right]\\ {A_{34}}\left( {{q_z}} \right) = {\partial _x}\left[ {M{\partial _z}{q_z}} \right] \end{array} \right. \end{eqnarray*} (A5) \begin{eqnarray*} \left\{ \begin{array}{@{}l@{}} {A_{41}}\left( {{v_x}} \right) = {\partial _z}\left( {\alpha M{\partial _x}{v_x}} \right)\\ {A_{42}}\left( {{v_z}} \right) = {\partial _z}\left( {\alpha M{\partial _z}{v_z}} \right)\\ {A_{43}}\left( {{q_x}} \right) = {\partial _z}\left( {M{\partial _x}{q_x}} \right)\\ {A_{44}}\left( {{q_z}} \right) = {\partial _z}\left( {M{\partial _z}{q_z}} \right) \end{array} \right.. \end{eqnarray*} (A6) APPENDIX B: SECOND-ORDER SPATIAL DERIVATIVES IN THE ROTATED CO-ORDINATE FRAME For the sake of simplicity, we let \begin{eqnarray*} \cos \theta = C,\,\,\sin \theta = S. \end{eqnarray*} (B1) Then, we have \begin{eqnarray*} \displaystyle\frac{\partial }{{\partial \bar{x}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial \bar{x}}}} \right] = \left\{ \begin{array}{@{}l@{}} CC\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] + CS\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right]\\ + CS\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] + SS\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right] \end{array} \right. \end{eqnarray*} (B2) \begin{eqnarray*} \displaystyle\frac{\partial }{{\partial \bar{z}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial \bar{z}}}} \right] = \left\{ \begin{array}{@{}l@{}} SS\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] - CS\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right]\\ - CS\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] + CC\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right] \end{array} \right. \end{eqnarray*} (B3) \begin{eqnarray*} \displaystyle\frac{\partial }{{\partial \bar{x}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial \bar{z}}}} \right] = \left\{ \begin{array}{@{}l@{}} - CS\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] + CC\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right]\\ - SS\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] + CS\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right] \end{array} \right. \end{eqnarray*} (B4) \begin{eqnarray*} \displaystyle\frac{\partial }{{\partial \bar{z}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial \bar{x}}}} \right] = \left\{ \begin{array}{@{}l@{}} - CS\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] - SS\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right]\\ + CC\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right] + CS\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left[ {f\left( {x,z} \right)\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right] \end{array} \right.. \end{eqnarray*} (B5) APPENDIX C: OPERATORS USED IN EQ. (18) Let us define the coefficient vectors |$T{R_k}$| (⁠|$k$| = 1,2,3,4) as \begin{eqnarray*} \begin{array}{@{}*{8}{l}@{}} {T{R_1} = \left[ {\begin{array}{*{20}{c}} {{C^2}}&{CS}&{CS}&{{S^2}} \end{array}} \right]}&{T{R_2} = \left[ {\begin{array}{*{20}{c}} { - CS}&{{C^2}}&{ - {S^2}}&{CS} \end{array}} \right]}\\ {T{R_3} = \left[ {\begin{array}{*{20}{c}} { - CS}&{ - {S^2}}&{{C^2}}&{CS} \end{array}} \right]}&{T{R_4} = \left[ {\begin{array}{*{20}{c}} {{S^2}}&{ - CS}&{ - CS}&{{C^2}} \end{array}} \right]} \end{array}. \end{eqnarray*} (C1) Here, |$C$| and |$S$| are defined in (B1). We also define the operator vectors |$DR( {m,v} )$| (⁠|$m$| represents material properties and |$v$| means wavefield) . \begin{eqnarray*} DR\left( {m,v} \right) = {\left[ {\begin{array}{*{20}{c}} {\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left( {m\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right),}&{\displaystyle\frac{\partial }{{\partial x{^{\prime}}}}\left( {m\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right),}&{\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left( {m\displaystyle\frac{{\partial v}}{{\partial x{^{\prime}}}}} \right),}&{\displaystyle\frac{\partial }{{\partial z{^{\prime}}}}\left( {m\displaystyle\frac{{\partial v}}{{\partial z{^{\prime}}}}} \right)} \end{array}} \right]^T}. \end{eqnarray*} (C2) The superscript ‘|$T$|’ means transpose operation. For example, |$T{R_4} \cdot DR( {\mu ,{v_x}} ) = SS\frac{\partial }{{\partial x{^{\prime}}}}[ {\mu \frac{{\partial {v_x}}}{{\partial x{^{\prime}}}}} ] - CS\frac{\partial }{{\partial x{^{\prime}}}}[ {\mu \frac{{\partial {v_x}}}{{\partial z{^{\prime}}}}} ] - CS\frac{\partial }{{\partial z{^{\prime}}}}[ {\mu \frac{{\partial {v_x}}}{{\partial x{^{\prime}}}}} ] + CC\frac{\partial }{{\partial z{^{\prime}}}}[ {\mu \frac{{\partial {v_x}}}{{\partial z{^{\prime}}}}} ]$|⁠. Then we get the concise expressions as \begin{eqnarray*} \left. \begin{array}{@{}l@{}} {{\bar{A}}_{11}}\left( {{v_x}} \right) = T{R_1} \cdot DR\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3},{v_x}} \right) + T{R_4} \cdot DR\left( {\mu ,{v_x}} \right)\\ {{\overset{\scriptscriptstyle\rightharpoonup}{A}}_{12}}\left( {{v_z}} \right) = T{R_2} \cdot DR\left( {{K^{sat}} - {{2\mu } \mathord{\left/ {\vphantom {{2\mu } 3}} \right. } 3},{v_z}} \right) + T{R_3} \cdot DR\left( {\mu ,{v_z}} \right)\\ {{\overset{\scriptscriptstyle\rightharpoonup}{A}}_{13}}\left( {{q_x}} \right) = T{R_1} \cdot DR\left( {\alpha M,{q_x}} \right)\\ {\bar{A}}_{14}\left( {{q_z}} \right) = T{R_2} \cdot DR\left( {\alpha M,{q_z}} \right) \end{array} \right\} \end{eqnarray*} (C3) \begin{eqnarray*} \left. \begin{array}{@{}l@{}} {{\bar{A}}_{21}}\left( {{v_x}} \right) = T{R_2} \cdot DR\left( {\mu ,{v_x}} \right) + TR{}_3 \cdot DR\left( {{K^{sat}} - {{2\mu } \mathord{\left/ {\vphantom {{2\mu } 3}} \right. } 3},{v_x}} \right)\\ {{\bar{A}}_{22}}\left( {{v_z}} \right) = T{R_1} \cdot DR\left( {\mu ,{v_z}} \right) + T{R_4} \cdot DR\left( {{K^{sat}} + {{4\mu } \mathord{\left/ {\vphantom {{4\mu } 3}} \right. } 3},{v_z}} \right)\\ {{\bar{A}}_{23}}\left( {{q_x}} \right) = TR{}_3 \cdot DR\left( {\alpha M,{q_x}} \right)\\ {{\bar{A}}_{24}}\left( {{q_z}} \right) = TR{}_4 \cdot DR\left( {\alpha M,{q_z}} \right) \end{array} \right\} \end{eqnarray*} (C4) \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{l}@{}} {{{\bar{A}}_{31}}\left( {{v_x}} \right) = TR{}_1 \cdot DR\left( {\alpha M,{v_x}} \right),}&{{{\bar{A}}_{32}}\left( {{v_z}} \right) = TR{}_2 \cdot DR\left( {\alpha M,{v_z}} \right)}\\ {{{\bar{A}}_{33}}\left( {{q_x}} \right) = TR{}_1 \cdot DR\left( {M,{q_x}} \right),}&{{{\bar{A}}_{34}}\left( {{q_z}} \right) = TR{}_2 \cdot DR\left( {M,{q_z}} \right)} \end{array}} \right\} \end{eqnarray*} (C5) \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{l}@{}} {{{\bar{A}}_{41}}\left( {{v_x}} \right) = = TR{}_3 \cdot DR\left( {\alpha M,{v_x}} \right),}&{{{\bar{A}}_{42}}\left( {{v_z}} \right) = TR{}_4 \cdot DR\left( {\alpha M,{v_z}} \right)}\\ {{{\bar{A}}_{43}}\left( {{q_x}} \right) = TR{}_3 \cdot DR\left( {M,{q_x}} \right),}&{{{\bar{A}}_{44}}\left( {{q_z}} \right) = TR{}_4 \cdot DR\left( {M,{q_z}} \right)} \end{array}} \right\}. \end{eqnarray*} (C6) APPENDIX D: THE |$R{^{\prime}}$| TERMS APPEARING IN EQ. (23) Let us define the operator vectors |${A{^{\prime}}_k}$| as \begin{eqnarray*} \left. {\begin{array}{@{}*{8}{c}@{}} {{{A{^{\prime}}}_1} = {{\left[ {\begin{array}{*{20}{c}} {{{A{^{\prime}}}_{11}}}&{{{A{^{\prime}}}_{12}}}&{{{A{^{\prime}}}_{21}}}&{{{A{^{\prime}}}_{22}}} \end{array}} \right]}^T}}&{{{A{^{\prime}}}_2} = {{\left[ {\begin{array}{*{20}{c}} {{{A{^{\prime}}}_{13}}}&{{{A{^{\prime}}}_{14}}}&{{{A{^{\prime}}}_{23}}}&{{{A{^{\prime}}}_{24}}} \end{array}} \right]}^T}}\\ {{{A{^{\prime}}}_3} = {{\left[ {\begin{array}{*{20}{c}} {{{A{^{\prime}}}_{31}}}&{{{A{^{\prime}}}_{32}}}&{{{A{^{\prime}}}_{41}}}&{{{A{^{\prime}}}_{41}}} \end{array}} \right]}^T}}&{{{A{^{\prime}}}_4} = {{\left[ {\begin{array}{*{20}{c}} {{{A{^{\prime}}}_{33}}}&{{{A{^{\prime}}}_{34}}}&{{{A{^{\prime}}}_{43}}}&{{{A{^{\prime}}}_{44}}} \end{array}} \right]}^T}} \end{array}} \right\}. \end{eqnarray*} (D1) Here |${A{^{\prime}}_{ij}}$| has the same form as |${A_{ij}}$| except the spatial derivatives are applied along the axis in the rotated coordinate system, e.g. |${A{^{\prime}}_{13}}( {{q_1}} ) = {\partial _{1{^{\prime}}}}[ {\alpha M{\partial _{1{^{\prime}}}}{q_1}} ]$|⁠. Then, with the definition of the coefficient vectors |$T{R_k}$|(⁠|$k$| = 1,2,3,4) (defined in C1), the rotated system terms in eq. (23) show even more concise expressions: \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{{R{^{\prime}}}_{11}}\left( {{v_x}} \right) = T{R_1} \cdot {{A{^{\prime}}}_1}\left( {{v_x}} \right)}&{{{R{^{\prime}}}_{12}}\left( {{v_z}} \right) = T{R_2} \cdot {{A{^{\prime}}}_1}\left( {{v_z}} \right)}\\ {{{R{^{\prime}}}_{13}}\left( {{q_x}} \right) = T{R_1} \cdot {{A{^{\prime}}}_1}\left( {{q_x}} \right)}&{{{R{^{\prime}}}_{14}}\left( {{q_z}} \right) = T{R_2} \cdot {{A{^{\prime}}}_2}\left( {{q_z}} \right)} \end{array}} \right\} \end{eqnarray*} (D2) \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{{R{^{\prime}}}_{21}}\left( {{v_x}} \right) = T{R_3} \cdot {{A{^{\prime}}}_1}\left( {{v_x}} \right)}&{{{R{^{\prime}}}_{22}}\left( {{v_z}} \right) = T{R_4} \cdot {{A{^{\prime}}}_1}\left( {{v_z}} \right)}\\ {{{R{^{\prime}}}_{23}}\left( {{q_x}} \right) = T{R_3} \cdot {{A{^{\prime}}}_2}\left( {{q_x}} \right)}&{{{R{^{\prime}}}_{24}}\left( {{q_z}} \right) = T{R_4} \cdot {{A{^{\prime}}}_2}\left( {{q_z}} \right)} \end{array}} \right\} \end{eqnarray*} (D3) \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{{R{^{\prime}}}_{31}}\left( {{v_x}} \right) = T{R_1} \cdot {{A{^{\prime}}}_3}\left( {{v_x}} \right)}&{{{R{^{\prime}}}_{32}}\left( {{v_z}} \right) = T{R_2} \cdot {{A{^{\prime}}}_3}\left( {{v_z}} \right)}\\ {{{R{^{\prime}}}_{33}}\left( {{q_x}} \right) = T{R_1} \cdot {{A{^{\prime}}}_4}\left( {{q_x}} \right)}&{{{R{^{\prime}}}_{34}}\left( {{q_z}} \right) = T{R_2} \cdot {{A{^{\prime}}}_4}\left( {{q_z}} \right)} \end{array}} \right\} \end{eqnarray*} (D4) \begin{eqnarray*} \left. {\begin{array}{@{}*{2}{c}@{}} {{{R{^{\prime}}}_{41}}\left( {{v_x}} \right) = T{R_3} \cdot {{A{^{\prime}}}_3}\left( {{v_x}} \right)}&{{{R{^{\prime}}}_{42}}\left( {{v_z}} \right) = T{R_4} \cdot {{A{^{\prime}}}_3}\left( {{v_z}} \right)}\\ {{{R{^{\prime}}}_{43}}\left( {{q_x}} \right) = T{R_3} \cdot {{A{^{\prime}}}_4}\left( {{q_x}} \right)}&{{{R{^{\prime}}}_{44}}\left( {{q_z}} \right) = T{R_4} \cdot {{A{^{\prime}}}_4}\left( {{q_z}} \right)} \end{array}} \right\}. \end{eqnarray*} (D5) APPENDIX E: PORE PRESSURE CALCULATION The Biot constitutive equation for pore pressure is given in the frequency-domain by \begin{eqnarray*} {p_f} = \frac{1}{{i\omega }}\left( {\alpha M\nabla \bullet {{\bf v}} + M\nabla \bullet {{\bf q}}} \right). \end{eqnarray*} (E1) Here |$\nabla \bullet {{\bf v}} = {\partial _x}{v_x} + {\partial _z}{v_z}$|and |$\nabla \bullet {{\bf q}} = {\partial _x}{q_x} + {\partial _z}{q_z}$|⁠. In the classic coordinate system \begin{eqnarray*} \begin{array}{*{20}{c}} {{\partial _x}{v_x} = \frac{1}{{2\Delta x}}\left[ {{v_x}\left( {m + 1,n} \right) - {v_x}\left( {m - 1,n} \right)} \right]}&{{\partial _z}{v_z} = \frac{1}{{2\Delta z}}\left[ {{v_z}\left( {m,n + 1} \right) - {v_z}\left( {m,n - 1} \right)} \right]} \end{array} \end{eqnarray*} (E2) \begin{eqnarray*} {\partial _x}{v_x} + {\partial _z}{v_z} = \frac{1}{{2\Delta x}}\left[ {{v_x}\left( {m + 1,n} \right) - {v_x}\left( {m - 1,n} \right)} \right] + \frac{1}{{2\Delta z}}\left[ {{v_z}\left( {m,n + 1} \right) - {v_z}\left( {m,n - 1} \right)} \right]. \end{eqnarray*} (E3) In the rotated coordinate system, using eq. (16) we have \begin{eqnarray*} {\partial _{\bar{x}}}{v_x} = \frac{{\partial {v_x}}}{{\partial x{^{\prime}}}}\frac{{\partial x{^{\prime}}}}{{\partial x}} + \frac{{\partial {v_x}}}{{\partial z{^{\prime}}}}\frac{{\partial z{^{\prime}}}}{{\partial x}} = \cos \theta \frac{{\partial {v_x}}}{{\partial x{^{\prime}}}} + \sin \theta \frac{{\partial {v_x}}}{{\partial z{^{\prime}}}} \end{eqnarray*} (E4) \begin{eqnarray*} {\partial _{\bar{z}}}{v_z} = \frac{{\partial {v_z}}}{{\partial x{^{\prime}}}}\frac{{\partial x{^{\prime}}}}{{\partial z}} + \frac{{\partial {v_z}}}{{\partial z{^{\prime}}}}\frac{{\partial z{^{\prime}}}}{{\partial z}} = - \sin \theta \frac{{\partial {v_z}}}{{\partial x{^{\prime}}}} + \cos \theta \frac{{\partial {v_z}}}{{\partial z{^{\prime}}}} \end{eqnarray*} (E5) \begin{eqnarray*} \frac{{\partial {v_x}}}{{\partial x{^{\prime}}}} = \frac{1}{{2\sqrt {{{\left( {\Delta x} \right)}^2} + {{\left( {\Delta z} \right)}^2}} }}\left[ {{v_x}\left( {m + 1,n - 1} \right) - {v_x}\left( {m - 1,n + 1} \right)} \right] \end{eqnarray*} (E6) \begin{eqnarray*} \frac{{\partial {v_x}}}{{\partial z{^{\prime}}}} = \frac{1}{{2\sqrt {{{\left( {\Delta x} \right)}^2} + {{\left( {\Delta z} \right)}^2}} }}\left[ {{v_x}\left( {m + 1,n + 1} \right) - {v_x}\left( {m - 1,n - 1} \right)} \right] \end{eqnarray*} (E7) \begin{eqnarray*} {\partial _{\bar{x}}}{v_x} + {\partial _{\bar{z}}}{v_z} = \cos \theta \frac{{\partial {v_x}}}{{\partial x{^{\prime}}}} + \sin \theta \frac{{\partial {v_x}}}{{\partial z{^{\prime}}}} - \sin \theta \frac{{\partial {v_z}}}{{\partial x{^{\prime}}}} + \cos \theta \frac{{\partial {v_z}}}{{\partial z{^{\prime}}}} \end{eqnarray*} (E8) \begin{eqnarray*} {\partial _{\bar{x}}}{v_x} + {\partial _{\bar{z}}}{v_z} = \frac{1}{{2\sqrt {{{\left( {\Delta x} \right)}^2} + {{\left( {\Delta z} \right)}^2}} }}\left\{ {\begin{array}{@{}*{1}{r}@{}} {\cos \theta \left[ {{v_x}\left( {m + 1,n - 1} \right) - {v_x}\left( {m - 1,n + 1} \right)} \right]}\\ { + \sin \theta \left[ {{v_x}\left( {m + 1,n + 1} \right) - {v_x}\left( {m - 1,n - 1} \right)} \right]}\\ { - \sin \theta \left[ {{v_z}\left( {m + 1,n - 1} \right) - {v_z}\left( {m - 1,n + 1} \right)} \right]}\\ { + \cos \theta \left[ {{v_z}\left( {m + 1,n + 1} \right) - {v_z}\left( {m - 1,n - 1} \right)} \right]} \end{array}} \right.. \end{eqnarray*} (E9) Here |$\cos \theta$| and |$\sin \theta$| given in eq. (16). By applying eqs (E3) and (E9), the pore pressure is numerically calculated with the mix-grid approach as \begin{eqnarray*} \left. {\begin{array}{@{}*{1}{r}@{}} {{p_f}\left( {m,n} \right) = \displaystyle\frac{a}{{i\omega }}\left[ {\alpha \left( {m,n} \right)M\left( {m,n} \right)\left( {{\partial _x}{v_x} + {\partial _z}{v_z}} \right) + M\left( {m,n} \right)\left( {{\partial _x}{q_x} + {\partial _z}{q_z}} \right)} \right]}\\ { + \displaystyle\frac{{1 - a}}{{i\omega }}\left[ {\alpha \left( {m,n} \right)M\left( {m,n} \right)\left( {{\partial _{\bar{x}}}{v_x} + {\partial _{\bar{z}}}{v_z}} \right) + M\left( {m,n} \right)\left( {{\partial _{\bar{x}}}{q_x} + {\partial _{\bar{z}}}{q_z}} \right)} \right]} \end{array}} \right\}. \end{eqnarray*} (E10) Similar to the eq. (26), |$a$| is set to be 0.5461 (Jo et al. 1996; Hustedt et al. 2004). © The Author(s) 2018. Published by Oxford University Press on behalf of The Royal Astronomical Society. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Frequency-domain seismic wave modelling in heterogeneous porous media using the mixed-grid finite-difference method
JF - Geophysical Journal International
DO - 10.1093/gji/ggy410
DA - 2019-01-01
UR - https://www.deepdyve.com/lp/oxford-university-press/frequency-domain-seismic-wave-modelling-in-heterogeneous-porous-media-jo6uveEFVV
SP - 34
VL - 216
IS - 1
DP - DeepDyve
ER -