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

Authors: Jinghuai Gao; Yijie Zhang
Journal: Mathematical Problems in Engineering
 (2013)

---

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


---

## Body

## 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 Section2, 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
(1)∂vx∂t=R22(∂τxx∂x+∂τxz∂z)-R12∂s∂x+B2(wx-vx),(2)∂vz∂t=R22(∂τxz∂x+∂τzz∂z)-R12∂s∂z+B2(wz-vz),(3)∂wx∂t=R11∂s∂x-R12(∂τxx∂x+∂τxz∂z)-B1(wx-vx),(4)∂wz∂t=R11∂s∂z-R12(∂τxz∂x+∂τzz∂z)-B1(wz-vz),(5)∂τxx∂t=P∂vx∂x+(P-2N)∂vz∂z+Q(∂wx∂x+∂wz∂z)+qx(wx-vx)+qz(wz-vz),(6)∂τzz∂t=(P-2N)∂vx∂x+P∂vz∂z+Q(∂wx∂x+∂wz∂z)+qx(wx-vx)+qz(wz-vz),(7)∂τxz∂t=N(∂vx∂z+∂vz∂x),(8)∂s∂t=Q(∂vx∂x+∂vz∂z)+R(∂wx∂x+∂wz∂z)+rx(wx-vx)+rz(wz-vz),
where (vx,vz) is the solid particle velocity vector, (wx,wz) is the fluid particle velocity vector, and (τxx,τxz,τzz) is the solid stress tensor. s is related to the fluid pressure p and the porosity ϕ according to
(9)s=-ϕp.
The coefficients P, Q, R, and N are Biot’s elastic constants. P and N correspond to the familiar Lame coefficients in the theory of elasticity. The coefficient R 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 Q 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 coefficientb is related to Darcy’s coefficient of permeability κ by
(10)b=μϕ2κ,
where μ is the fluid viscosity.The other parameters in (1)–(8) are defined as Rij=ρij/(ρ11ρ22-ρ122), B1=b(R11+R12), B2=b(R12+R22), qx=(Q/ϕ)(∂ϕ/∂x), qz=(Q/ϕ)(∂ϕ/∂z), rx=(R/ϕ)(∂ϕ/∂x), and rz=(R/ϕ)(∂ϕ/∂z).In staggered-grid scheme, particle velocity components, solid stress components, and fluid pressure are located at different points, which may be discretized as(vx)p,qn, (wx)p,qn, (vz)p+1/2,q+1/2n, (wz)p+1/2,q+1/2n, (τxx)p+1/2,qn-1/2, (τzz)p+1/2,qn-1/2, (s)p+1/2,qn-1/2, and (τxz)p,q+1/2n-1/2. The subscripts refer to the spatial indices, and the superscripts refer to the time indices. Thus, with a grid spacing h and time step τ, for example, the expression (vx)p,qn represents the x component of the solid particle velocity evaluated at the point [x+ph,z+qh] and at the time t+nτ.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
(11)(vx)p,qn=Aei[kx(x+ph)+kz(z+qh)-ω(t+nτ)],(wx)p,qn=Cei[kx(x+ph)+kz(z+qh)-ω(t+nτ)],(vz)p+1/2,q+1/2n=Dei[kx(x+(p+1/2)h)+kz(z+(q+1/2)h)-ω(t+nτ)],(wz)p+1/2,q+1/2n=Eei[kx(x+(p+1/2)h)+kz(z+(q+1/2)h)-ω(t+nτ)],(τxx)p+1/2,qn+1/2=Fei[kx(x+(p+1/2)h)+kz(z+qh)-ω(t+(n+1/2)τ)],(τzz)p+1/2,qn+1/2=Gei[kx(x+(p+1/2)h)+kz(z+qh)-ω(t+(n+1/2)τ)],(s)p+1/2,qn+1/2=Hei[kx(x+(p+1/2)h)+kz(z+qh)-ω(t+(n+1/2)τ)],(τxz)p,q+1/2n+1/2=Kei[kx(x+ph)+kz(z+(q+1/2)h)-ω(t+(n+1/2)τ)],
where ω is the angular frequency, (kx,kz) is the wavenumber vector, and i=-1.Substituting(A.1)–(A.20) of the Appendix and (11) into (1)–(8) and simplifying,
(12)gA=R22(fxF+fzK)-R12fxH,gD=R22(fxK+fzG)-R12fzH,gC=R11fxH-R12(fxF+fzK),gE=R11fzH-R12(fxK+fzG),gF=PfxA+Q(fxC+fzE)+(P-2N)fzD,gG=(P-2N)fxA+Q(fxC+fzE)+PfzD,gK=N(fxD+fzA),gH=Q(fxA+fzD)+R(fxC+fzE)
with
(13)g=-sin(ωτ2),fx=τh∑m=1Mamsin((m-12)kxh),fz=τh∑m=1Mamsin((m-12)kzh),
where am(m=1,2,…,M) are the SFD coefficients [28] and 2M corresponds to the order of accuracy.Equations (12) can be rewritten into matrix form(14)[0000R22fx0R22fz-R12fx00000R22fzR22fx-R12fz0000-R12fx0-R12fzR11fx00000-R12fz-R12fxR11fzPfx(P-2N)fzQfxQfz0000(P-2N)fxPfzQfxQfz0000NfzNfx000000QfxQfzRfxRfz0000][ADCEFGKH]=g[ADCEFGKH].That is, Wξ=gξ, g is the eigenvalue of the matrix W, we obtain
(15)g2=0,(16)g2=R22N(fx2+fz2),(17)g2=(Z+W2)(fx2+fz2),(18)g2=(Z-W2)(fx2+fz2),
where Z=(1/2)R11R+(1/2)R22P-R12Q and W=R2R112-2PRR11R22-4RQR11R12+P2R222-4PQR22R12+4PRR122+4Q2R11R22.For porous media, Biot [1] derived the velocities of fast P wave, S wave, and slow P wave. The velocity of S wave is
(19)vS=[Nρ11(1-ρ122/(ρ11ρ22))]1/2.
In order to derive the velocities of fast P wave and slow P wave, we introduce a reference velocity vc defined by
(20)vc2=Hρ.
The following parameters are introduced:
(21)σ11=PH,σ22=RH,σ12=QH,γ11=ρ11ρ,γ22=ρ22ρ,γ12=ρ12ρ,
where H=P+R+2Q and ρ=ρ11+2ρ12+ρ22.For an equation ofz(22)(σ11σ22-σ122)z2-(σ11γ22+σ22γ11-2σ12γ12)z+(γ11γ22-γ122)=0,
there are two roots z1 and z2 corresponding to two velocities v1,v2,
(23)v12=vc2z1,(24)v22=vc2z2.
The larger velocity between v1 and v2 is the velocity of fast P wave vP_fast, and the lower is the velocity of slow P wave vP_slow.Comparing (16)–(18) with (19), (23), and (24), (16)–(18) can be expressed as
(25)g2=vS2(fx2+fz2),(26)g2=vPslow2(fx2+fz2),(27)g2=vPfast2(fx2+fz2),
where vP_fast, vP_slow, and vS 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 kx=kcosθ and kz=ksinθ, where θ is a propagation direction angle of the plane wave. Dispersion relation of the minimum velocity can be written as
(28)(∑m=1Mamsin((m-0.5)khcosθ))2+(∑m=1Mamsin((m-0.5)khsinθ))2=(r-1sin(0.5ωτ))2,
where r=vτ/h and v=min(vS,vP_slow,vP_fast).Grid dispersion parameter in SFD modeling is defined by using (28),
(29)δwave(θ)=vFDwavevwave=2rwavekhsin-1(rwaveq)
with
(30)q=(∑m=1Mamsin((m-0.5)khcosθ))2+(∑m=1Mamsin((m-0.5)khsinθ))2,
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.Figure1(a) shows the variation of the dispersion parameter δ with kh 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 kh is related to velocity, we can change kh into f, which is independent of velocity. We plot the variation of the dispersion parameter δ with f in Figure 1(b). It can be concluded that the numerical dispersion for high-frequency components reduces with the increase of velocity for fixed M and h. As kh ranges from 0 to π (at Nyquist frequency), f changes from 0 to v/2h. 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.Dispersion curves for the same order of accuracy but different velocities. (a) Variation of dispersion parameterδ with kh, (b) variation of dispersion parameter δ with f, (c) variation of dispersion parameter log10|ε| with f. Here, v=500, 1500, 2500, 3500, 4500 m/s, τ=0.0001 s, h=2 m, and M=10. When log10|ε|<-16, we set log10|ε|=-16.
(a)
(b)
(c)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,(31)ε=hvFD-hv=hv(δ-1-1).
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,M, and f. 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 frequencyfmax and the maximum error η, the following inequality is satisfied:
(32)|ε(v,M,f)|≤η,whenf≤fmax.
The maximum frequency fmax is related to the source frequency, and the maximum error η should be larger than the maximum value of ε.Figure2 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, h=2 m, fmax=30 Hz, and η=10-9. The values of M 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,v=500, 1500, 2500, 3500, 4500 m/s, τ=0.0001 s, h=2 m, fmax=30 Hz, η=10-9, and M=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 Figure3. 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 h, and the maximum frequency fmax.Figure 3
Orders of accuracy for different velocities and maximum errors. Here,τ=0.0001 s, h=2 m, fmax=30 Hz. η=10-10, 10-9, 10-8, 10-7, 10-6, and 10-5. The values of M 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 Table1. We use a 30 Hz Ricker wavelet in time domain and Gaussian function in space domain located at x=400 m, z=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.
ρ
12
ρ
11
ρ
22
P
Q
R
N
b
ϕ
v
P_fast
v
S
v
P_slow
(kg/m3)
(N/m2)
(Nsm4)
(m/s)
1
-
260
1947
520
5.3
e
9
5.2
e
8
3.8
e
8
1.8
e
9
1
e
3
0.25
1856
1005
733
2
-
104
2174
208
7.6
e
9
2.4
e
8
1.4
e
8
2.6
e
9
1
e
3
0.1
1950
1100
764
3
-
104
2384
156
9.2
e
9
2.1
e
8
3.4
e
7
3.1
e
9
1
e
3
0.05
2030
1162
514
4
-
31
2446
62
1.2
e
10
1.1
e
8
3.9
e
7
4.1
e
9
1
e
3
0.03
2276
780
1298
5
-
330
2092
660
7.8
e
9
9.5
e
8
5.3
e
8
2.9
e
9
1
e
3
0.33
2215
1232
718
6
-
100
2430
255
5.5
e
10
7
e
8
7.9
e
8
1.2
e
10
7
e
5
0.25
4860
2240
1740Figure 4
The minimum velocities for the horizontally layered porous model and the source location. A30 Hz Ricker wavelet in time domain and Gaussian function in space domain located at x=400 m, z=360 m is used to generate the vibration. Model size is 800m×800m, and h=2 m, τ=0.0001 s.The orders of accuracy for this model are presented in Table2. They are determined by (32) with η=10-9, fmax=30 Hz, τ=0.0001 s, and h=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.
Minimum velocity (m/s)
514
718
733
764
780
1740
M
6
4
4
4
4
3The snapshot ofz 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 z components of the solid particle velocity calculated by the SFD method with fixed- and variable-order accuracy, respectively.Snapshots of thez 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 (M=6) and the SFD method with variable-order accuracy (M=6,3,3,3,3,2). No absorbing boundary conditions are used in the modeling.
(a)
The conventional method
(b)
The new method
(c)
The conventional method
(d)
The new methodThe modeling records of thez components of the solid particle velocity at x=400 m, z=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) Modeling records and (b) the difference of modeling records of thez components of the solid particle velocity at x=400 m, z=400 m for the horizontally layered porous model by SFD schemes with fixed- and variable-order accuracy.
(a)
(b)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 Table3. 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.
Type of SFDoperators
Value ofM
CPU time of 1000 recursion times (s)
Fixed
6
688
Variable
6, 4, 4, 4, 4, 3
572Data in Table4 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 Table1. The difference of these models is the thickness of Layer 3 and Layer 6. Values of M=6,4,4,4,4,3 are used for six different velocities.
Thickness (m)
CPU time of 1000 recursion times (s)
Layer 3
Layer 6
Fixed
Variable
280
240
692
606
40
240
688
572
40
440
698
558The 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 withkh or f 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.



---

*Source: 157071-2013-05-09.xml*