Staggered-Grid Finite Difference Method with Variable-Order Accuracy for Porous Media

The numerical modeling of wave field in porous media generally requires more computation time than that of acoustic or elastic media. Usually used finite difference methods adopt finite difference operators with fixed-order accuracy to calculate space derivatives for a heterogeneous medium. A finite difference scheme with variable-order accuracy for acoustic wave equation has been proposed to reduce the computation time. In this paper, we develop this scheme for wave equations in porous media based on dispersion relation with high-order staggered-grid finite difference (SFD) method. High-order finite difference operators are adopted for low-velocity regions, and low-order finite difference operators are adopted for high-velocity regions.Dispersion analysis and modeling results demonstrate that the proposed SFD method can decrease computational costs without reducing accuracy.


Introduction
Seismic wave modeling for porous media is usually used to study the properties of rocks and to characterize the seismic response of geologic formation.Furthermore, the theory about poroelasticity is a suitable characterization for the environment of hydrocarbon reservoirs.
The most popular theory about poroelasticity is developed by Biot [1][2][3], which is the basis for solving the wave propagation problems in porous media.A variety of different numerical methods have been used for poroelasticity modeling [4,5], such as spectral method [6], finite difference method [7], time domain method [8], discontinuous Galerkin method [9], and finite volume method [10].Owing to the low memory requirement and computational cost, finite difference methods [11,12] are widely applied for porous wave equations.
Compared with the numerical modeling of acoustic or elastic wave propagation, the porous wave propagation may spend more computation time.As we know, the wavelength in a lower velocity medium is smaller than that in a highvelocity medium with a fixed source frequency [24] and the accuracy of high-velocity regions is higher than that of low-velocity regions with a fixed-order finite difference operator for heterogeneous media.In addition, high-order finite difference operators provide better accuracy than loworder finite difference operators.Therefore, we can adopt high-order finite difference operators for low-velocity regions and low-order finite difference operators for high-velocity regions to reduce the computational costs.
The remainder of this paper is organized as follows.In Section 2, Biot's equations are reformulated into the firstorder, velocity-stress system; we derive the dispersion relation based on high-order staggered-grid finite difference method for porous media.In Section 3, we propose a method to determine the orders of accuracy for different velocities automatically.The validity of this method is demonstrated by dispersion analysis.In Section 4, we use the staggeredgrid finite difference method with fixed-and variable-order accuracy to simulate porous media.The modeling results demonstrate the efficiency of our method.We state the conclusion of this paper in Section 5.

Dispersion Relation for Porous Media
Biot established the dynamic equations in a porous elastic solid saturated by a compressible viscous fluid, and predicted that two dilatational waves and one rotational wave exist in the fluid-saturated porous solid.A first-order hyperbolic system [25][26][27] is developed, which is equivalent to Biot's equations, whose vector of unknowns consists of the solid and fluid particle velocity components, the solid stress components, and the fluid pressure.According to Biot's theory, the equations of motion for 2D porous media are given by where (V  , V  ) is the solid particle velocity vector, (  ,   ) is the fluid particle velocity vector, and (  ,   ,   ) is the solid stress tensor. is related to the fluid pressure  and the porosity  according to The coefficients , , , and  are Biot's elastic constants. and  correspond to the familiar Lame coefficients in the theory of elasticity.The coefficient  is a measure of the pressure required on the fluid to force a certain volume of the fluid into the aggregate, while the total volume remains constant.The coefficient  is the nature of a coupling between the volume change of solid and that of the fluid.And  11 ,  12 , and  22 are mass coefficients which take into account the fact that the relative fluid flow through the pores is not uniform.
The coefficient  is related to Darcy's coefficient of permeability  by where  is the fluid viscosity.In staggered-grid scheme, particle velocity components, solid stress components, and fluid pressure are located at different points, which may be discretized as +1/2, , and (  ) −1/2 ,+1/2 .The subscripts refer to the spatial indices, and the superscripts refer to the time indices.Thus, with a grid spacing ℎ and time step , for example, the expression (V  )  , represents the  component of the solid particle velocity evaluated at the point [ + ℎ,  + ℎ] and at the time  + .
The Appendix gives the discretization of the derivatives in (1)-( 8) based on high-order staggered-grid finite difference scheme.Using the plane wave theory, we let where  is the angular frequency, (  ,   ) is the wavenumber vector, and  = √ −1.11) into (1)-( 8) and simplifying,

Substituting (A.1)-(A.20) of the Appendix and (
where   ( = 1, 2, . . ., ) are the SFD coefficients [28] and 2 corresponds to the order of accuracy.Equations ( 12) can be rewritten into matrix form That is, W = ,  is the eigenvalue of the matrix W, we obtain where  = (1/2) For porous media, Biot [1] derived the velocities of fast P wave, S wave, and slow P wave.The velocity of S wave is In order to derive the velocities of fast P wave and slow P wave, we introduce a reference velocity V  defined by The following parameters are introduced: there are two roots  1 and  2 corresponding to two velocities The larger velocity between V 1 and V 2 is the velocity of fast P wave V P fast , and the lower is the velocity of slow P wave V P slow .

Mathematical Problems in Engineering
Comparing ( 16)-( 18) with ( 19), (23), and ( 24), ( 16)-( 18) can be expressed as where V P fast , V P slow , and V S are the velocities of fast P wave, slow P wave, and S wave, respectively.Equations ( 25)-( 27) are the dispersion relations based on high-order staggered-grid finite difference method for porous media.

SFD Method with Variable-Order Accuracy
The source frequency is usually fixed in the modeling so that the wavelength in a low-velocity region is smaller than that in a high-velocity region.For a heterogeneous model, the accuracy of high-velocity regions is higher than that of lowvelocity regions with a fixed-order finite difference operator to calculate the space derivatives.Moreover, high-order finite difference operator provides better accuracy than that of low-order finite difference operator.Therefore, we can adopt high-order finite difference operators for low-velocity regions and low-order finite difference operators for high-velocity regions.Note that three velocities exist in each grid for porous media, and we should choose the minimum velocity among them to determine the order of accuracy.Substitute ( 13) into ( 25), (26), or (27), and let   =  cos  and   =  sin , where  is a propagation direction angle of the plane wave.Dispersion relation of the minimum velocity can be written as where  = V/ℎ and V = min(V S , V P slow , V P fast ).Grid dispersion parameter in SFD modeling is defined by using (28), where the subscript wave represents slow P wave, S wave, or fast P wave.If  in (29) equals unity, there is no dispersion.If  is larger or less than unity, the dispersion will occur.
Figure 1(a) shows the variation of the dispersion parameter  with ℎ for the same order of accuracy but different velocities.From this figure it can be noted that the dispersion curves are nearly consistent for different velocities.As the independent variable ℎ is related to velocity, we can change ℎ into , which is independent of velocity.We plot the variation of the dispersion parameter  with  in Figure 1(b).It can be concluded that the numerical dispersion for highfrequency components reduces with the increase of velocity for fixed  and ℎ.As ℎ ranges from 0 to  (at Nyquist frequency),  changes from 0 to V/2ℎ.For different velocities, the values of maximum frequency are different.Note that the source frequency in the modeling is predetermined and is usually correlated with the minimum velocity for fixed-grid size.Therefore, we can reduce the order of accuracy for the high-velocity regions and achieve nearly the same accuracy as the high-order finite difference operators for the low-velocity regions.
The difference between SFD propagation time and exact propagation time through one grid could describe the error due to the numerical dispersion in SFD scheme, If  equals zero, there is no dispersion; when  is much greater or less than zero, large dispersion will occur.Following the definition of SFD error,  is related to V, , and . Figure 1(c) illustrates that the SFD error for low velocity is larger than that for high velocity with a certain frequency, which implies that we can reduce the order of accuracy for the high-velocity regions without reducing the accuracy.
Here, we adopt a method [24] to determine the orders of accuracy for different velocities adaptively.
For the given maximum frequency  max and the maximum error , the following inequality is satisfied: The maximum frequency  max is related to the source frequency, and the maximum error  should be larger than the maximum value of . Figure 2 shows dispersion curves for different orders of accuracy and different velocities; the involved parameters are V = 500, 1500, 2500, 3500, 4500 m/s,  = 0.0001 s, ℎ = 2 m,  max = 30 Hz, and  = 10 −9 .The values of  determined by (32) are 6, 3, 3, 2, and 2 for the five different velocities, respectively.From this figure it can be noted that the loworder finite difference operators can be adopted for highvelocity regions to obtain an accuracy that is better than or equal to the accuracy of the high-order finite difference operators for low-velocity regions.
The orders of accuracy with different velocities and maximum errors are presented in Figure 3.The figure demonstrates that the orders of accuracy generally decrease with the increase of velocities.In addition, the orders of accuracy determined by (32) are dependent on the time step , the grid spacing ℎ, and the maximum frequency  max .

Modeling Results
We adopt the high-order staggered-grid finite difference method for a horizontally layered porous model, whose parameters are shown in Table 1.We use a 30 Hz Ricker wavelet in time domain and Gaussian function in space domain located at  = 400 m,  = 360 m to generate the vibration, and the model size is 800 m × 800 m (see Figure 4).The orders of accuracy for this model are presented in Table 2.They are determined by (32) with  = 10 −9 ,  max = 30 Hz,  = 0.0001 s, and ℎ = 2 m.Note that highorder finite difference operators provide better accuracy than low-order finite difference operators.Therefore, the highest order of accuracy used in the variable-order accuracy SFD method is adopted as the order of accuracy for the fixedorder accuracy SFD method which can be served as reference solution.
The snapshot of  component of the fluid particle velocity calculated by the SFD method with variable-order accuracy (see Figure 5(b)) is almost the same as Figure 5(a), which is calculated by the SFD method with fixed-order accuracy.The same conclusion can be obtained from the comparison of Figures 5(c) and 5(d), which are the snapshots of  components of the solid particle velocity calculated by the SFD method with fixed-and variable-order accuracy, respectively.
The modeling records of the  components of the solid particle velocity at  = 400 m,  = 400 m calculated by the SFD method with fixed-and variable-order accuracy are presented in Figure 6(a).The difference between them is shown in Figure 6(b).From this figure, it can be noted that the modeling results calculated by the two methods are almost identical, which demonstrate the validity of the staggeredgrid finite difference method with variable-order accuracy.
Computational efficiency of the SFD method with variable-order accuracy for this model is demonstrated by CPU time (HP with an Intel Q8400 Core 2 quard CPU and 4.00 GB of memory), shown in Table 3.It can be seen that the SFD method with variable-order accuracy can save about 17% of the CPU time in comparison with the SFD method with fixed-order accuracy.Data in Table 4 show the CPU time for different models.We can conclude that the efficiency depends on the characteristic of the models.That is, the savings of CPU time increase with the increase of high-velocity regions and the decrease of low-velocity regions.
The modeling results and the CPU time demonstrate that the SFD method with variable-order accuracy results in a decrease in computation time.

Conclusion
We have developed a staggered-grid finite difference scheme with variable-order accuracy for porous media.We use a method based on dispersion relation to determine the orders of accuracy for different velocity regions.The variations of dispersion parameters with ℎ or  imply the validity of the new scheme.The comparison of numerical modeling results and CPU time between the SFD schemes with fixed-and variable-order accuracy indicates that the proposed scheme can enhance the computational efficiency without reducing the accuracy.

Figure 4 :
Figure 4: The minimum velocities for the horizontally layered porous model and the source location.A 30 Hz Ricker wavelet in time domain and Gaussian function in space domain located at  = 400 m,  = 360 m is used to generate the vibration.Model size is 800 m × 800 m, and ℎ = 2 m,  = 0.0001 s.

Figure 5 :
Figure 5: Snapshots of the  components of fluid (a)-(b) and solid (c)-(d) particle velocity at 192 ms for the horizontally layered porous model by the SFD method with fixed-order accuracy ( = 6) and the SFD method with variable-order accuracy ( = 6, 3, 3, 3, 3, 2).No absorbing boundary conditions are used in the modeling.

Figure 6 :
Figure 6: (a) Modeling records and (b) the difference of modeling records of the  components of the solid particle velocity at  = 400 m,  = 400 m for the horizontally layered porous model by SFD schemes with fixed-and variable-order accuracy.

Table 1 :
Parameters for the horizontally layered porous model.

Table 2 :
Orders of accuracy for different velocities with the same maximum frequency.