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

Abstract

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.

1. 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–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.

To improve the efficiency and accuracy, several variants of finite difference methods have been investigated to simulate porous media [13, 14]. These include implicit finite difference method [15], variable-grids [16], irregular-grids, staggered-grids [17], discontinuous-grids [18], quadrangle-grids [19], rotated-staggered-grids [20] finite difference methods, and spatially varying time steps finite difference methods [21]. In order to avoid the spurious reflection from the artificial boundaries in finite difference modeling, perfectly matched layer (PML) absorbing boundary conditions [22, 23] have been widely used.

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 high-velocity 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 low-order 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 first-order, 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 staggered-grid 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.

2. 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–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 𝜕𝑣_𝑥𝜕𝑡=𝑅₂₂𝜕𝜏_𝑥𝑥+𝜕𝑥𝜕𝜏_𝑥𝑧𝜕𝑧−𝑅₁₂𝜕𝑠𝜕𝑥+𝐵₂𝑤_𝑥−𝑣_𝑥,𝜕𝑣_𝑧𝜕𝑡=𝑅₂₂𝜕𝜏_𝑥𝑧+𝜕𝑥𝜕𝜏_𝑧𝑧𝜕𝑧−𝑅₁₂𝜕𝑠𝜕𝑧+𝐵₂𝑤_𝑧−𝑣_𝑧,𝜕𝑤_𝑥𝜕𝑡=𝑅₁₁𝜕𝑠𝜕𝑥−𝑅₁₂𝜕𝜏_𝑥𝑥+𝜕𝑥𝜕𝜏_𝑥𝑧𝜕𝑧−𝐵₁𝑤_𝑥−𝑣_𝑥,𝜕𝑤_𝑧𝜕𝑡=𝑅₁₁𝜕𝑠𝜕𝑧−𝑅₁₂𝜕𝜏_𝑥𝑧+𝜕𝑥𝜕𝜏_𝑧𝑧𝜕𝑧−𝐵₁𝑤_𝑧−𝑣_𝑧,𝜕𝜏_𝑥𝑥𝜕𝑡=𝑃𝜕𝑣_𝑥+()𝜕𝑥𝑃−2𝑁𝜕𝑣_𝑧𝜕𝑧+𝑄𝜕𝑤_𝑥+𝜕𝑥𝜕𝑤_𝑧𝜕𝑧+𝑞_𝑥𝑤_𝑥−𝑣_𝑥+𝑞_𝑧𝑤_𝑧−𝑣_𝑧,𝜕𝜏_𝑧𝑧=()𝜕𝑡𝑃−2𝑁𝜕𝑣_𝑥𝜕𝑥+𝑃𝜕𝑣_𝑧𝜕𝑧+𝑄𝜕𝑤_𝑥+𝜕𝑥𝜕𝑤_𝑧𝜕𝑧+𝑞_𝑥𝑤_𝑥−𝑣_𝑥+𝑞_𝑧𝑤_𝑧−𝑣_𝑧,𝜕𝜏_𝑥𝑧𝜕𝑡=𝑁𝜕𝑣_𝑥+𝜕𝑧𝜕𝑣_𝑧,𝜕𝑥𝜕𝑠𝜕𝑡=𝑄𝜕𝑣_𝑥+𝜕𝑥𝜕𝑣_𝑧𝜕𝑧+𝑅𝜕𝑤_𝑥+𝜕𝑥𝜕𝑤_𝑧𝜕𝑧+𝑟_𝑥𝑤_𝑥−𝑣_𝑥+𝑟_𝑧𝑤_𝑧−𝑣_𝑧, (1)(2)(3)(4)(5)(6)(7)(8) where (𝑣_𝑥,𝑣_𝑧) 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 𝑠=−𝜙𝑝. (9) 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 𝜌₁₁, 𝜌₁₂, and 𝜌₂₂ 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 𝑏=𝜇𝜙²𝜅, (10) where 𝜇 is the fluid viscosity.

The other parameters in (1)–(8) are defined as 𝑅_𝑖𝑗=𝜌_𝑖𝑗/(𝜌₁₁𝜌₂₂−𝜌²₁₂), 𝐵₁=𝑏(𝑅₁₁+𝑅₁₂), 𝐵₂=𝑏(𝑅₁₂+𝑅₂₂), 𝑞_𝑥=(𝑄/𝜙)(𝜕𝜙/𝜕𝑥), 𝑞_𝑧=(𝑄/𝜙)(𝜕𝜙/𝜕𝑧), 𝑟_𝑥=(𝑅/𝜙)(𝜕𝜙/𝜕𝑥), and 𝑟_𝑧=(𝑅/𝜙)(𝜕𝜙/𝜕𝑧).

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,𝑞+1/2}, (𝑤_𝑧)^𝑛_{𝑝+1/2,𝑞+1/2}, (𝜏_𝑥𝑥)^𝑛−1/2_{𝑝+1/2,𝑞}, (𝜏_𝑧𝑧)^𝑛−1/2_{𝑝+1/2,𝑞}, (𝑠)^𝑛−1/2_{𝑝+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 (𝑣_𝑥)^𝑛_𝑝,𝑞 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 𝑣_𝑥^𝑛_𝑝,𝑞=𝐴𝑒^{𝑖[𝑘𝑥(𝑥+𝑝ℎ)+𝑘𝑧(𝑧+𝑞ℎ)−𝜔(𝑡+𝑛𝜏)]},𝑤_𝑥^𝑛_𝑝,𝑞=𝐶𝑒^{𝑖[𝑘𝑥(𝑥+𝑝ℎ)+𝑘𝑧(𝑧+𝑞ℎ)−𝜔(𝑡+𝑛𝜏)]},𝑣_𝑧^𝑛_{𝑝+1/2,𝑞+1/2}=𝐷𝑒^{𝑖[𝑘𝑥(𝑥+(𝑝+1/2)ℎ)+𝑘𝑧(𝑧+(𝑞+1/2)ℎ)−𝜔(𝑡+𝑛𝜏)]},𝑤_𝑧^𝑛_{𝑝+1/2,𝑞+1/2}=𝐸𝑒^{𝑖[𝑘𝑥(𝑥+(𝑝+1/2)ℎ)+𝑘𝑧(𝑧+(𝑞+1/2)ℎ)−𝜔(𝑡+𝑛𝜏)]},𝜏_𝑥𝑥^𝑛+1/2_{𝑝+1/2,𝑞}=𝐹𝑒^{𝑖[𝑘𝑥(𝑥+(𝑝+1/2)ℎ)+𝑘𝑧(𝑧+𝑞ℎ)−𝜔(𝑡+(𝑛+1/2)𝜏)]},𝜏_𝑧𝑧^𝑛+1/2_{𝑝+1/2,𝑞}=𝐺𝑒^{𝑖[𝑘𝑥(𝑥+(𝑝+1/2)ℎ)+𝑘𝑧(𝑧+𝑞ℎ)−𝜔(𝑡+(𝑛+1/2)𝜏)]},(𝑠)^𝑛+1/2_{𝑝+1/2,𝑞}=𝐻𝑒^{𝑖[𝑘𝑥(𝑥+(𝑝+1/2)ℎ)+𝑘𝑧(𝑧+𝑞ℎ)−𝜔(𝑡+(𝑛+1/2)𝜏)]},𝜏_𝑥𝑧^𝑛+1/2_{𝑝,𝑞+1/2}=𝐾𝑒^{𝑖[𝑘𝑥(𝑥+𝑝ℎ)+𝑘𝑧(𝑧+(𝑞+1/2)ℎ)−𝜔(𝑡+(𝑛+1/2)𝜏)]}, (11) where 𝜔 is the angular frequency, (𝑘_𝑥,𝑘_𝑧) is the wavenumber vector, and 𝑖=√−1.

Substituting (A.1)–(A.20) of the Appendix and (11) into (1)–(8) and simplifying, 𝑔𝐴=𝑅₂₂𝑓_𝑥𝐹+𝑓_𝑧𝐾−𝑅₁₂𝑓_𝑥𝐻,𝑔𝐷=𝑅₂₂𝑓_𝑥𝐾+𝑓_𝑧𝐺−𝑅₁₂𝑓_𝑧𝐻,𝑔𝐶=𝑅₁₁𝑓_𝑥𝐻−𝑅₁₂𝑓_𝑥𝐹+𝑓_𝑧𝐾,𝑔𝐸=𝑅₁₁𝑓_𝑧𝐻−𝑅₁₂𝑓_𝑥𝐾+𝑓_𝑧𝐺,𝑔𝐹=𝑃𝑓_𝑥𝐴+𝑄𝑓_𝑥𝐶+𝑓_𝑧𝐸+(𝑃−2𝑁)𝑓_𝑧()𝑓𝐷,𝑔𝐺=𝑃−2𝑁_𝑥𝐴+𝑄𝑓_𝑥𝐶+𝑓_𝑧𝐸+𝑃𝑓_𝑧𝑓𝐷,𝑔𝐾=𝑁_𝑥𝐷+𝑓_𝑧𝐴,𝑔𝐻=𝑄𝑓_𝑥𝐴+𝑓_𝑧𝐷+𝑅𝑓_𝑥𝐶+𝑓_𝑧𝐸 (12) with 𝑔=−sin𝜔𝜏2,𝑓_𝑥=𝜏ℎ^𝑀_𝑚=1𝑎_𝑚1sin𝑚−2𝑘_𝑥ℎ,𝑓_𝑧=𝜏ℎ^𝑀_𝑚=1𝑎_𝑚1sin𝑚−2𝑘_𝑧ℎ, (13) where 𝑎_𝑚(𝑚=1,2,…,𝑀) are the SFD coefficients [28] and 2𝑀 corresponds to the order of accuracy.

Equations (12) can be rewritten into matrix form0000𝑅₂₂𝑓_𝑥0𝑅₂₂𝑓_𝑧−𝑅₁₂𝑓_𝑥00000𝑅₂₂𝑓_𝑧𝑅₂₂𝑓_𝑥−𝑅₁₂𝑓_𝑧0000−𝑅₁₂𝑓_𝑥0−𝑅₁₂𝑓_𝑧𝑅₁₁𝑓_𝑥00000−𝑅₁₂𝑓_𝑧−𝑅₁₂𝑓_𝑥𝑅₁₁𝑓_𝑧𝑃𝑓_𝑥(𝑃−2𝑁)𝑓_𝑧𝑄𝑓_𝑥𝑄𝑓_𝑧()𝑓0000𝑃−2𝑁_𝑥𝑃𝑓_𝑧𝑄𝑓_𝑥𝑄𝑓_𝑧0000𝑁𝑓_𝑧𝑁𝑓_𝑥000000𝑄𝑓_𝑥𝑄𝑓_𝑧𝑅𝑓_𝑥𝑅𝑓_𝑧𝐴𝐷𝐶𝐸𝐹𝐺𝐾𝐻𝐴𝐷𝐶𝐸𝐹𝐺𝐾0000𝐻=𝑔. (14)That is, 𝐖𝝃=𝑔𝝃, 𝑔 is the eigenvalue of the matrix 𝐖, we obtain 𝑔²𝑔=0,₂=𝑅₂₂𝑁𝑓²_𝑥+𝑓²_𝑧𝑔₂,=√𝑍+𝑊2𝑓²_𝑥+𝑓²_𝑧,𝑔²=√𝑍−𝑊2𝑓²_𝑥+𝑓²_𝑧, (15)(16)(17)(18) where 𝑍=(1/2)𝑅₁₁𝑅+(1/2)𝑅₂₂𝑃−𝑅₁₂𝑄 and 𝑊=𝑅²𝑅²₁₁−2𝑃𝑅𝑅₁₁𝑅₂₂−4𝑅𝑄𝑅₁₁𝑅₁₂+𝑃²𝑅²₂₂−4𝑃𝑄𝑅₂₂𝑅₁₂+4𝑃𝑅𝑅²₁₂+4𝑄²𝑅₁₁𝑅₂₂.

For porous media, Biot [1] derived the velocities of fast P wave, S wave, and slow P wave. The velocity of S wave is 𝑣_S=𝑁𝜌₁₁(1−𝜌²₁₂/(𝜌₁₁𝜌₂₂))^1/2. (19) In order to derive the velocities of fast P wave and slow P wave, we introduce a reference velocity 𝑣_𝑐 defined by 𝑣²_𝑐=𝐻𝜌. (20) The following parameters are introduced: 𝜎₁₁=𝑃𝐻,𝜎₂₂=𝑅𝐻,𝜎₁₂=𝑄𝐻,𝛾₁₁=𝜌₁₁𝜌,𝛾₂₂=𝜌₂₂𝜌,𝛾₁₂=𝜌₁₂𝜌, (21) where 𝐻=𝑃+𝑅+2𝑄 and 𝜌=𝜌₁₁+2𝜌₁₂+𝜌₂₂.

For an equation of 𝑧𝜎₁₁𝜎₂₂−𝜎²₁₂𝑧²−𝜎₁₁𝛾₂₂+𝜎₂₂𝛾₁₁−2𝜎₁₂𝛾₁₂𝑧+𝛾₁₁𝛾₂₂−𝛾²₁₂=0, (22) there are two roots 𝑧₁ and 𝑧₂ corresponding to two velocities 𝑣₁,𝑣₂, 𝑣²₁=𝑣²_𝑐𝑧₁,𝑣²₂=𝑣²_𝑐𝑧₂. (23)(24) The larger velocity between 𝑣₁ and 𝑣₂ is the velocity of fast P wave 𝑣_{P_fast}, and the lower is the velocity of slow P wave 𝑣_{P_slow}.

Comparing (16)–(18) with (19), (23), and (24), (16)–(18) can be expressed as 𝑔²=𝑣²_S𝑓²_𝑥+𝑓²_𝑧𝑔₂,=𝑣²_Pslow𝑓²_𝑥+𝑓²_𝑧𝑔₂,=𝑣²_Pfast𝑓²_𝑥+𝑓²_𝑧, (25)(26)(27) where 𝑣_{P_fast}, 𝑣_{P_slow}, and 𝑣_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.

3. 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 low-velocity 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 ^𝑀_𝑚=1𝑎_𝑚sin((𝑚−0.5)𝑘ℎcos𝜃)²+^𝑀_𝑚=1𝑎_𝑚sin((𝑚−0.5)𝑘ℎsin𝜃)²=𝑟⁻¹sin(0.5𝜔𝜏)², (28) where 𝑟=𝑣𝜏/ℎ and 𝑣=min(𝑣_S,𝑣_{P_slow},𝑣_{P_fast}).

Grid dispersion parameter in SFD modeling is defined by using (28), 𝛿_wave(𝜃)=𝑣_FDwave𝑣_wave=2𝑟_wave𝑘ℎsin⁻¹𝑟_wave√𝑞 (29) with 𝑞=^𝑀_𝑚=1𝑎_𝑚sin((𝑚−0.5)𝑘ℎcos𝜃)²+^𝑀_𝑚=1𝑎_𝑚sin((𝑚−0.5)𝑘ℎsin𝜃)², (30) 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 high-frequency components reduces with the increase of velocity for fixed 𝑀 and ℎ. As 𝑘ℎ ranges from 0 to 𝜋 (at Nyquist frequency), 𝑓 changes from 0 to 𝑣/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.

(a)

(b)

(c)

Figure 1 

Dispersion curves for the same order of accuracy but different velocities. (a) Variation of dispersion parameter 𝛿 with 𝑘ℎ, (b) variation of dispersion parameter 𝛿 with 𝑓, (c) variation of dispersion parameter log10|𝜀| with 𝑓. Here, 𝑣=500, 1500, 2500, 3500, 4500 m/s, 𝜏=0.0001 s, ℎ=2 m, and 𝑀=10. When log10|𝜀|<−16, we set log10|𝜀|=−16.

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, ℎ𝜀=𝑣_FD−ℎ𝑣=ℎ𝑣𝛿⁻¹−1. (31) 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 𝑣,𝑀, 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: 𝜀𝑣,𝑀,𝑓≤𝜂,when𝑓≤𝑓_max. (32) 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 𝑣=500, 1500, 2500, 3500, 4500 m/s, 𝜏=0.0001 s, ℎ=2 m, 𝑓_max=30 Hz, and 𝜂=10⁻⁹. 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 low-order finite difference operators can be adopted for high-velocity 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.

Figure 2 

Dispersion curves for different orders of accuracy and velocities. Here, 𝑣=500, 1500, 2500, 3500, 4500 m/s, 𝜏=0.0001 s, ℎ=2 m, 𝑓max=30 Hz, 𝜂=10−9, and 𝑀=6, 3, 3, 2, and 2 determined by (32). When log10|𝜀|<−16, we set log10|𝜀|=−16.

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.

Figure 3 

Orders of accuracy for different velocities and maximum errors. Here, 𝜏=0.0001 s, ℎ=2 m, 𝑓max=30 Hz. 𝜂=10−10, 10−9, 10−8, 10−7, 10−6, and 10−5. The values of 𝑀 are determined by (32). When log10|𝜀|<−16, we set log10|𝜀|=−16.

4. 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 800m×800m (see Figure 4).

Table 1 

Parameters for the horizontally layered porous model.

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 800m×800m, and ℎ=2 m, 𝜏=0.0001 s.

The orders of accuracy for this model are presented in Table 2. They are determined by (32) with 𝜂=10⁻⁹, 𝑓_max=30 Hz, 𝜏=0.0001 s, and ℎ=2 m. Note that high-order 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 fixed-order accuracy SFD method which can be served as reference solution.

Table 2 

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

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.

(a) The conventional method

(b) The new method

(c) The conventional method

(d) The new method

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.

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 staggered-grid finite difference method with variable-order accuracy.

(a)

(b)

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.

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.

Table 3 

CPU time of SFD modeling for the horizontally layered model.

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.

Table 4 

CPU time of SFD modeling for different models, whose parameters are shown in Table 1. The difference of these models is the thickness of Layer 3 and Layer 6. Values of 𝑀=6,4,4,4,4,3 are used for six different velocities.

The modeling results and the CPU time demonstrate that the SFD method with variable-order accuracy results in a decrease in computation time.

5. 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.

Appendix

In this appendix, we discretize the derivatives of wave equations based on high-order staggered-grid finite difference scheme for porous media.

The derivatives in (1) and (3) can be discretized as follows: 𝜕𝑣_𝑥=𝑣𝜕𝑡^{𝑥𝑛+1𝑝,𝑞}−𝑣_𝑥^𝑛_𝑝,𝑞𝜏,𝜕𝑤_𝑥=𝑤𝜕𝑡^{𝑥𝑛+1𝑝,𝑞}−𝑤_𝑥^𝑛_𝑝,𝑞𝜏,𝜕𝜏_𝑥𝑥=1𝜕𝑥ℎ^𝑀_𝑚=1𝑎_𝑚𝜏_𝑥𝑥()^𝑛+1/2_{𝑝+𝑚−1/2,𝑞}−𝜏_𝑥𝑥()^𝑛+1/2_{𝑝−𝑚−1/2,𝑞},𝜕𝜏_𝑥𝑧=1𝜕𝑧ℎ^𝑀_𝑚=1𝑎_𝑚𝜏_𝑥𝑧()^𝑛+1/2_{𝑝,𝑞+𝑚−1/2}−𝜏_𝑥𝑧()^𝑛+1/2_{𝑝,𝑞−𝑚−1/2},𝜕𝑠=1𝜕𝑥ℎ^𝑀_𝑚=1𝑎_𝑚(𝑠)₍₎^𝑛+1/2_{𝑝+𝑚−1/2,𝑞}−(𝑠)₍₎^𝑛+1/2_{𝑝−𝑚−1/2,𝑞}. (A.1)(A.2)(A.3)(A.4)(A.5)

Equations (2) and (4) can be discretized using the following equations: 𝜕𝑣_𝑧=𝑣𝜕𝑡^{𝑧𝑛+1𝑝+1/2,𝑞+1/2}−𝑣_𝑧^𝑛_{𝑝+1/2,𝑞+1/2}𝜏,𝜕𝑤_𝑧=𝑤𝜕𝑡^{𝑧𝑛+1𝑝+1/2,𝑞+1/2}−𝑤_𝑧^𝑛_{𝑝+1/2,𝑞+1/2}𝜏,𝜕𝜏_𝑥𝑧=1𝜕𝑥ℎ^𝑀_𝑚=1𝑎_𝑚𝜏_𝑥𝑧^𝑛+1/2_{𝑝+𝑚,𝑞+1/2}−𝜏_𝑥𝑧^𝑛+1/2_{𝑝−(𝑚−1),𝑞+1/2},𝜕𝜏_𝑧𝑧=1𝜕𝑧ℎ^𝑀_𝑚=1𝑎_𝑚𝜏_𝑧𝑧^𝑛+1/2_{𝑝+1/2,𝑞+𝑚}−𝜏_𝑧𝑧^𝑛+1/2_{𝑝+1/2,𝑞−(𝑚−1)},𝜕𝑠=1𝜕𝑧ℎ^𝑀_𝑚=1𝑎_𝑚(𝑠)^𝑛+1/2_{𝑝+1/2,𝑞+𝑚}−(𝑠)^𝑛+1/2_{𝑝+1/2,𝑞−(𝑚−1)}. (A.6)(A.7)(A.8)(A.9)(A.10)

The derivatives in (7) can be discretized as follows: 𝜕𝜏_𝑥𝑧=𝜏𝜕𝑡^{𝑥𝑧𝑛+1/2𝑝,𝑞+1/2}−𝜏_𝑥𝑧^𝑛−1/2_{𝑝,𝑞+1/2}𝜏,𝜕𝑣_𝑥=1𝜕𝑧ℎ^𝑀_𝑚=1𝑎_𝑚𝑣_𝑥^𝑛_{𝑝,𝑞+𝑚}−𝑣_𝑥^𝑛_{𝑝,𝑞−(𝑚−1)},𝜕𝑣_𝑧=1𝜕𝑥ℎ^𝑀_𝑚=1𝑎_𝑚𝑣_𝑧^𝑛_{𝑝+(𝑚−1/2),𝑞+1/2}−𝑣_𝑧^𝑛_{𝑝−(𝑚−1/2),𝑞+1/2}. (A.11)(A.12)(A.13)

Equations (5), (6), and (8) can be discretized using the following equations: 𝜕𝜏_𝑥𝑥=𝜏𝜕𝑡^{𝑥𝑥𝑛+1/2𝑝+1/2,𝑞}−𝜏_𝑥𝑥^𝑛−1/2_{𝑝+1/2,𝑞}𝜏,𝜕𝜏_𝑧𝑧=𝜏𝜕𝑡^{𝑧𝑧𝑛+1/2𝑝+1/2,𝑞}−𝜏_𝑧𝑧^𝑛−1/2_{𝑝+1/2,𝑞}𝜏,𝜕𝑠=(𝑠)𝜕𝑡^{𝑛+1/2𝑝+1/2,𝑞}−(𝑠)^𝑛−1/2_{𝑝+1/2,𝑞}𝜏,𝜕𝑣_𝑥=1𝜕𝑥ℎ^𝑀_𝑚=1𝑎_𝑚𝑣_𝑥^𝑛_{𝑝+𝑚,𝑞}−𝑣_𝑥^𝑛_{𝑝−(𝑚−1),𝑞},𝜕𝑤_𝑥=1𝜕𝑥ℎ^𝑀_𝑚=1𝑎_𝑚𝑤_𝑥^𝑛_{𝑝+𝑚,𝑞}−𝑤_𝑥^𝑛_{𝑝−(𝑚−1),𝑞},𝜕𝑣_𝑧=1𝜕𝑧ℎ^𝑀_𝑚=1𝑎_𝑚𝑣_𝑧^𝑛_{𝑝+1/2,𝑞+(𝑚−1/2)}−𝑣_𝑧^𝑛_{𝑝+1/2,𝑞−(𝑚−1/2)},𝜕𝑤_𝑧=1𝜕𝑧ℎ^𝑀_𝑚=1𝑎_𝑚𝑤_𝑧^𝑛_{𝑝+1/2,𝑞+(𝑚−1/2)}−𝑤_𝑧^𝑛_{𝑝+1/2,𝑞−(𝑚−1/2)}, (A.14)(A.15)(A.16)(A.17)(A.18)(A.19)(A.20) where ℎ is the grid size, 𝜏 is the time step, 𝑎_𝑚(1,2,…,𝑀) are staggered-grid finite difference coefficients.

Acknowledgment

This work is supported by the National Science and Technology Major Project of China (2011ZX05023-005-009, 2011ZX05044).