Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 157071, 10 pages

http://dx.doi.org/10.1155/2013/157071

Research Article

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

Jinghuai Gao1,2 and Yijie Zhang1,2

¹Institute of Wave and Information, Xi'an Jiaotong University, Xi'an 710049, China

²National Engineering Laboratory for Offshore Oil Exploration, Xi'an 710049, China

Received 9 January 2013; Accepted 7 April 2013

Academic Editor: Alex Elias-Zuniga

Copyright © 2013 Jinghuai Gao and Yijie Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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

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

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 , , , , , , , and . 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

where is the angular frequency, is the wavenumber vector, and .

Substituting – of the Appendix and (11) into (1)–(8) and simplifying,

with

where are the SFD coefficients [28] and corresponds to the order of accuracy.

Equations (12) can be rewritten into matrix form

That is, , is the eigenvalue of the matrix , we obtain

where and .

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 defined by

The following parameters are introduced:

where and .

For an equation of

there are two roots and corresponding to two velocities ,

The larger velocity between and is the velocity of fast P wave , and the lower is the velocity of slow P wave .

Comparing (16)–(18) with (19), (23), and (24), (16)–(18) can be expressed as

where , , and 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 and , where is a propagation direction angle of the plane wave. Dispersion relation of the minimum velocity can be written as

where and .

Grid dispersion parameter in SFD modeling is defined by using (28),

with

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 to (at Nyquist frequency), changes from to . 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 with . Here, , , , ,  m/s,  s,  m, and . When , we set .

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 , 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 and the maximum error , the following inequality is satisfied:

The maximum frequency 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 , , , ,  m/s,  s,  m,  Hz, and . The values of determined by (32) are , , , , and 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, , , , ,  m/s,  s,  m,  Hz, , and , , , , and determined by (32). When , we set .

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 .

Figure 3: Orders of accuracy for different velocities and maximum errors. Here,  s,  m,  Hz. , , , , , and . The values of are determined by (32). When , we set .

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  Hz Ricker wavelet in time domain and Gaussian function in space domain located at  m,  m to generate the vibration, and the model size is (see Figure 4).

Table 1: Parameters for the horizontally layered porous model.

(kg/m³)
(N/m²)
(Nsm⁴)

(m/s)

1

1947
520

0.25
1856
1005
733

2

2174
208

0.1
1950
1100
764

3

2384
156

0.05
2030
1162
514

4

2446
62

0.03
2276
780
1298

5

2092
660

0.33
2215
1232
718

6

2430
255

0.25
4860
2240
1740

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

The orders of accuracy for this model are presented in Table 2. They are determined by (32) with ,  Hz,  s, and  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

6
4
4
4
4
3

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  ms for the horizontally layered porous model by the SFD method with fixed-order accuracy () and the SFD method with variable-order accuracy (). No absorbing boundary conditions are used in the modeling.

The modeling records of the components of the solid particle velocity at  m,  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  m,  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.

Type of SFD

operators
Value of
CPU time of 1000 recursion times (s)

Fixed
6
688

Variable
6, 4, 4, 4, 4, 3
572

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 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
558

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:

Equations (2) and (4) can be discretized using the following equations:

The derivatives in (7) can be discretized as follows:

Equations (5), (6), and (8) can be discretized using the following equations:

where is the grid size, is the time step, 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).

References

M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low-frequency range,” Acoustical Society of America, vol. 28, pp. 168–178, 1956.

M. A. Biot, “Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Higher frequency range,” Acoustical Society of America, vol. 28, pp. 179–191, 1956.

M. A. Biot, “Mechanics of deformation and acoustic propagation in porous media,” vol. 33, pp. 1482–1498, 1962.

J. M. Carcione, Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic and Porous Media, Elsevier, Pergamon, 2001.

J. M. Carcione, C. Morency, and J. E. Santos, “Computational poroelasticity—a review,” Geophysics, vol. 75, no. 5, pp. 75A229–75A243, 2010.

G. Degrande and G. De Roeck, “FFT-based spectral analysis methodology for one-dimensional wave propagation in poroelastic media,” Transport in Porous Media, vol. 9, no. 1-2, pp. 85–97, 1992.

Y. J. Masson and S. R. Pride, “Finite-difference modeling of Biot's poroelastic equations across all frequencies,” Geophysics, vol. 75, no. 2, pp. N33–N41, 2010.

G. Chiavassa and B. Lombard, “Time domain numerical modeling of wave propagation in 2D heterogeneous porous media,” Journal of Computational Physics, vol. 230, no. 13, pp. 5288–5309, 2011.

B. Dupuy, L. De Barros, S. Garambois, and J. Virieux, “Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous Galerkin method,” Geophysics, vol. 76, no. 4, pp. N13–N28, 2011.

G. I. Lemoine, M. Y. Ou, and R. J. LeVeque, “High-resolution finite volume modeling of wave propagation in orthotropic poroelastic media,” SIAM Journal on Scientific Computing, vol. 35, no. 1, pp. B176–B206, 2013.

Z. Alterman and F. Karal, “Propagation of elastic waves in layered media by finite difference methods,” Bulletin of the Seismological Society of America, vol. 58, pp. 367–398, 1968.

R. M. Alford, K. R. Kelly, and D. M. Boore, “Accuracy of finite-difference modeling of the acoustic wave equation,” Geophysics, vol. 39, no. 6, pp. 834–841, 1974.

M. Kindelan, A. Kamel, and P. Sguazzero, “On the construction and efficiency of staggered numerical differentiators for the wave equation,” Geophysics, vol. 55, pp. 107–110, 1990.

O. Holberg, “Computational aspects of the choice of operator and sampling interval for numerical differentiation in large-scale simulation of wave phenomena,” Geophysical Prospecting, vol. 35, no. 6, pp. 629–655, 1987.

R. M. Beam and R. F. Warming, “An implicit finite-difference algorithm for hyperbolic systems in conservation-law form,” Journal of Computational Physics, vol. 22, no. 1, pp. 87–110, 1976.

P. S. Jensen, “Finite difference techniques for variable grids,” Computers and Structures, vol. 2, no. 1-2, pp. 17–29, 1972.

T. Özdenvar and G. A. McMechan, “Algorithms for staggered-grid computations for poroelastic, elastic, acoustic, and scalar wave equations,” Geophysical Prospecting, vol. 45, no. 3, pp. 403–420, 1997.

S. Aoi and H. Fujiwara, “3D finite-difference method using discontinuous grids,” Bulletin of the Seismological Society of America, vol. 89, no. 4, pp. 918–930, 1999.

J. Zhang, “Quadrangle-grid velocity-stress finite difference method for poroelastic wave equations,” Geophysical Journal International, vol. 139, no. 1, pp. 171–182, 1999.

G. S. O'Brien, “3D rotated and standard staggered finite-difference solutions to Biot's poroelastic wave equations: stability condition and dispersion analysis,” Geophysics, vol. 75, no. 4, pp. T111–T119, 2010.

E. Tessmer, “Seismic finite-difference modeling with spatially varying time steps,” Geophysics, vol. 65, no. 4, pp. 1290–1293, 2000.

Y. Q. Zeng, J. Q. He, and Q. H. Liu, “The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media,” Geophysics, vol. 66, no. 4, pp. 1258–1266, 2001.

S. D. Wang, “Perfectly matched layer absorbing boundary conditions for acoustic wave equation,” Oil Geophysical Prospecting, vol. 38, pp. 31–34, 2003 (Chinese).

Y. Liu and M. K. Sen, “Finite-difference modeling with adaptive variable-length spatial operators,” Geophysics, vol. 76, no. 4, pp. T79–T89, 2011.

J. Virieux, “P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method,” Geophysics, vol. 51, no. 4, pp. 889–901, 1986.

N. Dai, Finite difference simulation and imaging of seismic waves in complex media [Ph.D. thesis], University of Alberta, 1993.

N. Dai, A. Vafidis, and E. R. Kanasewich, “Wave propagation in heterogeneous, porous media: a velocity-stress, finite-difference method,” Geophysics, vol. 60, no. 2, pp. 327–340, 1995.

L. G. Dong, Z. T. Ma, J. Z. Cao et al., “A staggered-grid high-order difference method of one-order elastic wave equation,” Acta Geophysica Sinica, vol. 43, no. 3, pp. 411–419, 2000.