Finite-difference modeling with variable grid-size and adaptive time-step in porous media

1.   Introduction

Heterogeneous reservoirs, such as fractured or cavernous carbonate reservoirs, tight sandstone reservoirs, etc., have gradually become the hotspots of seismic exploration. These reservoirs exhibit two-phase features caused by the fluid in pores, cracks or fractures, and multi-scale features caused by complex structures and strong heterogeneous reservoirs. Both of these features should be considered in the explanation and prediction for heterogeneous oil/gas reservoirs. The numerical modeling of seismic waves in heterogeneous and porous reservoirs is an important tool for seismic interpretation (Wenzlau and Müller 2009).

Biot(1956, 1962) established the dynamic equations of elastic wave propagation in saturated porous media for the first time and predicted the existence of two types of longitudinal waves and one type of shear wave. The Biot theory has been regarded as a basis for wave propagation analysis in porous media. Since the analytical solutions of wave equations are difficult to obtain, different numerical methods are proposed by many authors. Finite-difference methods have been widely used because of low computational cost, ease of implementation, and strong adaptability. Staggered-grid finite difference (SGFD) method can further improve the precision and suppress dispersion without additional computational cost (Wu 2006). Hassanzadeh (1991) first solved the low-frequency Biot poroacoustic equations using a second-order finite-difference (FD) scheme. After that, FD method was extended and improved to solve poroelastic wave equations (Zhu and McMechan 1991; Dai et al. 1995; Carcione and Quiroga-Goode 1995; Carcione and Helle 1999; Santos et al. 2004; Wenzlau and Müller 2009; Sun and Yin 2011).

SGFD method uses regular grid to discretize seismic models and may lead to numerical dispersion until spatial grid-size is small enough. That means, for low velocity layers or heterogeneous bodies (such as multi-scale structures or interfaces with high contrast between physical parameters, etc.) surrounded by high-velocity homogeneous media, grid-size should be small enough to assure sufficient sampling with small time-step according to the FD stability condition. However, simulation with this small grid-size and time-step leads to oversampling in surrounding high-velocity homogeneous media and results in additional computational cost. To achieve the balance of spatial sampling in geological structures at different scales, Moczo (1989) first proposed the idea of varied grid-size to decrease computational cost. This idea have been improved for elastic wave simulation (Jastram and Tessmer 1994; Pitarka 1999; Hayashi et al. 2001; Wang et al. 2001), but the flexibility and accuracy of these methods are insufficient since the spatial grid-size variation is limited to single direction and grid-size ratio is restricted under three. The variable-grid FD schemes then have been developed with grid-size variation in both spatial dimensions (Zhu et al. 2007; Li et al. 2008; Guilbot 2008). However, the computational efficiency is still not good enough because of unchanged time step.

The local variable time-step methods were proposed to solve this problem (Falk et al. 1998; Tessmer 2000; Kang and Baag 2004; Huang and Dong 2009; Sun and Yin 2011). These methods solve the temporal oversampling problem with detailed description of local structures and refined temporal effectiveness. In this paper, SGFD method with variable spatial grid-size and variable time-step is extended to simulate elastic wave propagation in porous media. Small grid-size and time-step is used for heterogeneous bodies such as fractured or cavernous bodies, regions with high physical parameter contrasts, etc., while big grid-size and time-step are used for surrounding homogeneous media. Different FD coefficients and wavefield interpolation are used to realize the transition of wave propagation between regions discretized by different spatial grid-size.

2.   Methodology

According to Biot's theory (Biot 1956), combined with Newtonian dynamic equations and Darcy's law (Zhang, et al. 2004), the stress and velocity equations of wave propagation in fluid-saturated, statistically isotropic, nonuniform porous media can be derived as

where A, N, Q, R are elastic parameters of porous media; δ_ij is Kronecker function. σ_ij is solid stress, S is fluid effective strain; A_ij=ρ_ij(ρ₁₁ρ₂₂-ρ₁₂²)^-1, ρ_ij is mass parameters and are related to the mass parameters of solid and fluid per unit volume, ρ₁ and ρ₂: ρ₁₁+ρ₁₂=(1-ϕ)ρ₁, ρ₂₂+ρ₁₂=ϕρ₂; ϕ is porosity; b is dissipation parameter.

2.1   Staggered-grid finite-difference scheme based on PML boundary condition

Discretize the model into numbers of grids with certain grid-size, define model parameters and wavefield variables in the main grids and staggered grids which are staggered by half the grid spacing, as shown in Fig. 1.

Figure  1.  Definition of the main grids and staggered grids used in SGFD modeling

DownLoad: Full-Size Img PowerPoint

The FD scheme is also staggered in time. Let k, k±, k ± 1 represent time , respectively, substitute the FD formation of first-order time derivation into Eq. (1), we obtain the SGFD equations (take Eq. (1c)) for example and let i=j=₁) as

where L_x⁺, L_x^- denote forward and backward finite-difference operators along x direction, respectively.

The (2N)th-order SGFD operator derived by Taylor series expansion (Wu 2006) is given by

where C_i^(N)(i=₁, ₂, ⋯, N) are FD coefficients determined by the following equations:

where F_ij=(2j-1)^2i-1, i, j=1, 2, ⋯, N, E=(1, 0, ⋯, 0)^T. Combined with the definition of grids and model parameters, the forward and backward SGFD operators are

Truncating the unbounded media with artificial boundaries will introduce spurious reflections. The perfectly matched layer (PML) absorbs outgoing waves by imposing artificial absorptive layers outside the regular media (He 2008; Hastings et al. 1996). In this paper, the SGFD formula with PML was derived by modifying Eq. (1) with complex coordinates.

In frequency domain, introduce a complex factor by which the Cartesian coordinates are extended to complex coordinates (He 2008)

where n represents spatial coordinates; is attenuation factor independent of frequency, d_n > 0 inside PML and d_n = 0 outside PML. According to Eq. (6), the operator ∂/∂n can be expressed in terms of the complex coordinates

where a time dependent term n is implied. The PML formulation is to replace n in Eq. (1) by the corresponding complex coordinate Take the velocity component along x direction for example, we have

In order to simplify the PML equations, the variables are split along spatial coordinates

then Eq. (8) can be split into two equations

Transform the equations back to time–space domain and substitute the SGFD operators of first-order time and spatial derivatives into, we obtain

where (2+d_nΔt)^-1 and h_n=2Δt(2+d_nΔt)^-1. Similar equations can be obtained for other velocity and stress components. It is obvious that outside PML, the combination of Eqs. (9) and (11) degenerate into Eq. (2), which means the introduction of PML do not change the wavefields outside PML.

2.2   Staggered-grid finite-difference scheme with variable grid-size and time-step

Geological bodies at different scales are discretized with different grid-size. Small grid-size is used for small-scaled structures or fractured bodies and big grid-size is used for their surrounding media. There are two methods to deal with transition zone between small and big grid-size regions: interpolation method and variable coefficient method (Sun and Ding 2012). This paper introduces variable coefficient method which is of higher simulation precision.

Let Δx and Δz represent small grid-size along x and z direction, respectively, M Δx, M Δz (M∈N₊ is the ratio of big and small grid-size) represent big grid-size along x and z direction, respectively. High-order SGFD operators are used in both big and small grid-size regions. In transition zone (gray shaded area in Fig. 2), the FD coefficients should be recalculated to maintain the symmetry of SGFD operators and suppress artificial reflections arises by the interface of big and small grid-size region.

Figure  2.  Configuration of model discretion with spatial grids (gray shaded area represents the transition zone)

DownLoad: Full-Size Img PowerPoint

Take sixth-order SGFD operator with M = 3 for example, for first-order spatial derivations along x direction, grid nodes used to calculate the SGFD operators are shown in Fig. 3, the black stars represent the computing nodes (center of the SGFD operators), black circles represent nodes involved for the calculation. Formulations and coefficients of the FD operator of three nodes in the left of the interface (between big grid-size region and small grid-size region) remain unchanged, while the nodes involved for calculation should be symmetric about the center node. Formulations of the FD operator at two nodes in the right of the interface becomes

Figure  3.  Formulations of FD operators in transition zone

DownLoad: Full-Size Img PowerPoint

C_i^(N) is FD coefficients determined by

where , E=(1, 0, ⋯, 0)^T.

The stability condition for elastic-wave propagation holds also for the poroelastic case (Wenzlau and Müller 2009). The stability of the numerical scheme with (2N)th-order spatial and second-order time FD operator for elastic-wave propagation (Huang and Dong 2009) is

where V_{P, max} is the maximum P-wave velocity in big and small grid-size region. Equation (14) must be satisfied in big and small grid-size region, respectively. Variable spatial grid-size method would not make the stability condition more rigid (Huang and Dong 2009).

To satisfy the stability condition, time-step in small grid-size region is adaptively defined according to grid-size using Eq. (14). Let Δt and P Δt represent small and big time-step, respectively, first calculate wavefields using big grid-size M Δx, M Δz and time-step P Δt, save the field variables of the transition zone, calculate the wavefields in small grid-size region at this moment using linear interpolation method, then do time stepping with time-step Δt until time level P Δt is reached. The transition of wavefield form big time-step to small time-step is achieved. The time-step ratio P can be arbitrary positive integer adaptively determined by stability condition (14).

3.   Numerical examples

3.1   Homogeneous model of porous media

We first examine the accuracy of the numerical solutions. Compare the elastic wave modeling results of variable grid-size and time-step SGFD (VSGFD) method with the analytical solution proposed by Dai et al. (1995) in a 2-D homogeneous poroelastic model. Model parameters are the same as those of the uniform model of Dai et al. (1995) (Table 1) and the damping coefficient b is equal to zero. The size of model is 1000 × 1000 m². P-wave point source of Ricker time function with main frequency of 30 Hz locates at (500 m, 500 m). The big grid-size is 2 m, time-step is 0.1 ms; small grid-size is 0.4 m, arranges form 50 to 600 m along x direction and 200 to 400 m along z direction, small time-step is 0.02 ms.

Table  1.  Elastic parameters of homogeneous porous model

 | Show Table

DownLoad: CSV

The fast P-wave velocity is 3210 m/s, slow P-wave velocity is 1180 m/s, and the theoretical arrival time is 0.195 and 0.531 s, respectively. The comparison of numerical results and analytical solutions (both have been normalized) are in Figs. 4 and 5. The arrival time, amplitude and frequency of numerical solutions fit those of analytical solutions very well. The VSGFD method is accurate enough to simulate the wave propagation in poroelastic media and was further tested in the following sections.

Figure  4.  Z component wave record of solid phase [Receiver locates at (100 m, 20 m)]

DownLoad: Full-Size Img PowerPoint

Figure  5.  Z component wave record of fluid phase [Receiver locates at (100 m, 20 m)]

DownLoad: Full-Size Img PowerPoint

3.2   Cavernous model of porous media

The model contains two oil-bearing caves in the surrounding homogeneous media, the radius are 4 and 40 m, respectively. The parameters are listed in Table 2. The corresponding poroelastic parameters A, N, Q, R, ρ₁₁, ρ₂₂, and ρ₁₂ are calculated using rock physical theories with the method proposed by Geertsma and Smit (1961) and Berryman (1980).The big grid-size is 9 m, time-step is 1 ms; small grid-size is 3 m (inside dashed rectangle of Fig. 6). P-wave source with dominant frequency of 40 Hz locates at the middle of the model surface.

Table  2.  Parameters of porous model with caves

 | Show Table

DownLoad: CSV

Figure  6.  Porous model with caves

DownLoad: Full-Size Img PowerPoint

The modeling results of SGFD using big grid-size and time-step (BSGFD) could not simulate the reflected waves of the smaller cave (radius is 4 m), results of VSGFD is similar to SSGFD and the reflection of the smaller cave can be seen at around 0.25 ms. The computational time is 396.2, 1282.3, 7232.3 s by BSGFD, VSGFD, and SSGFD algorithm, respectively. The VSGFD method significantly improved the simulation efficiency while maintained the accuracy (Figs. 7, 8, 9).

Figure  7.  X component shot records of solid phase in porous model with caves

DownLoad: Full-Size Img PowerPoint

Figure  8.  Z component shot records of solid phase in porous model with caves (without direct wave)

DownLoad: Full-Size Img PowerPoint

Figure  9.  Z component wave motion of solid phase in porous model with caves [receiver locates at (800, 0)]

DownLoad: Full-Size Img PowerPoint

3.3   Low-velocity reservoir model of porous media

The model contains a low-velocity reservoir in surrounding high-velocity layer, the model parameters are listed in Table 3. The big grid-size is 10 m, time-step is 1 ms; small grid-size is 2 m (inside dashed rectangle of Fig. 10). Vertical force source with dominant frequency of 30 Hz locates at the middle of the model surface.

Table  3.  Parameters of low-velocity reservoir model

 | Show Table

DownLoad: CSV

Figure  10.  Porous model with caves

DownLoad: Full-Size Img PowerPoint

For the horizontal interface, reflection amplitude and arrival time simulated by BSGFD, VSGFD, SSGFD has little difference (Figs. 11, 12, 13, 14). Dispersion (can be seen after 0.65 s in Fig. 13) exists in the reflections of low-velocity reservoir simulated by BSGFD method. The SSGFD method can attenuate the dispersion but spent too much computational cost. The computational time is 55870.2 s while that spent by VSGFD method is 2991.6 s, that means the simulation efficiency is much improved. Simulation with small grid-size (2 m) oversamples the Surrounding high-velocity media and has higher requirements for computer performance while VSGFD method overcomes this limitation. The results of VSGFD method have negligible differences with results of SSGFD method.

Figure  11.  X component shot records of solid phase in porous model with low-velocity reservoir

DownLoad: Full-Size Img PowerPoint

Figure  12.  Z component shot records of solid phase in porous model with low-velocity reservoir

DownLoad: Full-Size Img PowerPoint

Figure  13.  Enlarged figures of shot records inside dashed rectangle of Fig. 12

DownLoad: Full-Size Img PowerPoint

Figure  14.  Z component wave motion of Figs. 14a and 12a–c [receiver locates at (1000, 0)]

DownLoad: Full-Size Img PowerPoint

Figure 15a shows the Z component shot records simulated for conventional solid single-phase media (using the velocities and densities in Table 3) using the same VSGFD method. Compared with Fig. 12b, we can found that, for the horizontal interface, reflected waves based on porous media theory and single-phase medium theory have reasonable concordance, this is because the first and second layer are designed to be ideal solid phase media (porosity and oil saturation is zero) and can be equivalent to single phase layers. Inside the reservoir, the relative movements between solid and fluid give rise to energy losses, result in attenuation with varying degrees of elastic waves with different frequencies.

Figure  15.  Z component shot records for single-phase media and the difference between Fig. 12a and Fig. 12b

DownLoad: Full-Size Img PowerPoint

4.   Discussion and conclusions

A velocity-stress, staggered-grid finite-difference scheme with variable grid-size and time-step has been implemented successfully to simulate wave propagation in heterogeneous, porous media. Variable spatial grid-size is chosen to insure sufficient sampling of structures at different scales or layers with high velocity contrasts, time-step is determined adaptively according to grid-size and highest velocity by stability condition of the finite-difference scheme. The proposed method can improve simulation efficiency while insure the accuracy and is more flexible to adapt to parameter changes in formation models.

Forward modeling based on poroelastic theory is more reasonable to simulate the influence of rock porosity and fluid properties on seismic wavefield, however there are still some issues need to be dealt with for practical application, including calculation of elastic parameters A, N, Q, R, ρ₁₁, ρ₂₂, and ρ₁₂, loading of source, treatment of different scales of heterogeneous bodies, the analysis of the complex wavefield, etc. This paper provides an idea to simulate elastic wave propagation and seismic responses while considering the properties of rock pores and fluid in heterogeneous porous media while intensive studies are still needed, especially in low-velocity oil/gas reservoirs with small structures.