Rock moduli estimation of inhomogeneous two-phase media with finite difference modeling algorithm

Digital rock physics (DRP) builds a bridge between pore-scale physical processes and the macroscopic physical properties of rock. Its key paradigm is to image and digitize the pore space and mineral matrix of rock, then to simulate the response of various physical fields and calculate the equivalent elastic parameters of rocks through digital rocks and mathematical methods. In this paper, a new approach is proposed to estimate the rock moduli of two-phase media with the staggered-grid high-order finite difference method (FDM) of the Biot theory based on DRP. This new method not only takes into consideration the impact of fluid on elastic wave propagation in two-phase media, but it is also easy to understand and implement, improving the calculation accuracy, requiring less memory and improving on the weaknesses of conventional rock physics experiments which are time consuming and expensive. In order to estimate the rock moduli, we establish two models, and the digital rock sample is embedded in one of those. Using this method, it is possible to model the dynamic wave propagation and measure the time delay of the peak amplitude caused by the inhomogeneous structure of the digital rock sample, with the receivers set at the bottom of the two models. The time delay allows us to estimate the effective velocity of both compressional and shear waves, and therefore calculate the rock moduli. Additionally, comparison between the numerical simulated results obtained through this method and experimental results indicates that they agree well. Comparison with the numerical simulated results obtained via another method tests and verifies the accuracy and feasibility of the new method. Also, the equivalence conditions between this new method and the various rock physics models are inferred.


Introduction
Rock from oil and gas reservoirs is a typical porous medium, the macroscopic physical properties of which are related to a variety of micro-factors (Liu et al 2013). Therefore, quantitatively characterizing the natural relations between the macroscopic physical properties and the microstructure is of significance to the exploration and development of oil and gas.
Reliable rock physics models have proved useful by reducing the range of possible interpretations. However, conventional rock experiments and theories have encountered difficulties in accurately predicting the properties of rocks, especially for unconventional reservoirs (Yin et al 2017). Rock physics experiments (Wyllie et al 1956(Wyllie et al , 1958 are inefficient, high cost and time consuming. It is difficult to study the effects of micro-factors on microstructure parameters and macroscopic properties of the rock using conventional rock experiments. Rock physics theories, which assume ideal rock structural features, are not consistent with real rock microstructures (Andrä et al 2013a(Andrä et al , 2013b. The empirical relationships are also concluded by a set of rocks. It is not suitable to apply these formulas to reservoir rocks in other areas (Mavko et al (1998)).
Digital rock physics (DRP) is the hot topic in rock physics (Liu et al 2013, Chen et al 2014. The aim of DRP is to discover, understand and model the relations between macroscopic physical properties and pore scale parameters (Andrä et al 2013a, 2013b, Zong et al 2012a, Yin et al 2015, Sun et al 2012. With the technological development of modern micro imaging, high performance computers and the numerical simulation algorithm, DRP technology is already gradually maturing (Madonna et al 2013). Using high resolution representations of the complex pore geometry, DRP reveals the real microstructure of the rocks. DRP can be used to not only estimate the effective elastic properties (EEP) of rock (e.g., elastic, transport, and electrical properties), but also to better understand the physical processes governing these properties (Saenger et al 2011).
Currently, digital cores have been extensively studied and applied in logging and oil-and-gas development, whereas there has been limited study and application in seismic exploration. The main numerical method used to estimate elastic moduli based on the 3D digital core is the finite element method (FEM). Arns et al (2002) calculated elastic parameters based on a Fontainebleau digital core with FEM, which corresponds to Gassmann (1951) theory. Grechka et al (2006a), Grechka and Kachanov (2006b), Grechka (2007) also used FEM to study the elastic properties of porous fractured media based on a digital core. Makarynska et al (2008) computationally simulated the EEP of partially fluidsaturated rocks, and verified that the numerical results are in agreement with those of the Gassmann-Wood (1955) equation at low frequency, and the Gassmann-Hill (1963) equation in high frequency if it is fluid-saturated rock. Sun et al (2015) applied FEM to study shale. Liu (2010) computationally simulated the effective elastic moduli with FEM based on a 3D digital core, and verified the reasonableness of the FEM. Zhao and Sun (2013) applied the FEM to analyze the effect of different methods of cementation on elastic properties and seepage characteristics of rock. However, the FEM has low computational efficiency (Yin et al 2016).
Another method used to estimate elastic moduli is the finite difference method (FDM). Saenger et al (2000) modeled the propagation of elastic waves using a new rotated staggered grid. For fluid-saturated porous media, Saenger et al (2004) calculated its equivalent EEP, which laid the foundation for DRP to be applied to seismic exploration. Saenger et al (2005) was concerned with numerical considerations of viscous fluid effects on wave propagation in porous media with FDM. Saenger et al (2007) used both the elastic and viscoelastic versions of the FDM algorithm to study gas hydrates and to simulate propagation of Biot's slow wave. Saenger et al (2011) studied the effect of pore fluid viscosity on EEP based on the digital rocks using the generalized Maxwell body and rotated staggered grid. Yin et al (2016) established the scalar wave equation and estimated porous media moduli with ordinary-grid FDM.
The numerical method in this paper is staggered-grid FDM, which defines velocity and stress in a staggered grid, respectively. The major difference between this and the ordinary grid is that the different components of velocity and stress field are not known at the same node (Madariaga 1976, Virieux (1986). The snapshots and simulated results of an actual model show that staggered-grid FDM is more accurate and can decrease the grid dispersion in conventional FDM. (Dong et al 2000, Pei andMu (2003)). Also, it is easy to understand and implement compared to rotated staggered grid FDM.
In this paper, a new approach is proposed to estimate the moduli of two-phase media by the staggered-grid high-order FDM based on DRP. This new method not only takes into consideration the impact of fluid on elastic wave propagation in two-phase media, it is also easy to understand and implement, improving the calculation accuracy, requiring less memory and improving on the weaknesses of the conventional rock physics experiment, which is time consuming and expensive. This paper establishes the first-order velocity and stress equations of 3D two-phase media based on Biot theory (Biot 1956). For the first-order velocity and stress equations of Biot theory, we use staggered-grid high-order FDM to solve the equations. The staggered-grid high-order FDM algorithm is the new approach to model dynamic wave propagation and obtain the effective velocity of both compressional and shear waves based on a 3D digital core. As a result, we have estimated the effective bulk moduli and shear moduli. Additionally, the equivalent conditions between this method and the theoretical rock physics models are inferred, and the effectiveness and rationality of this method has also been verified by comparing the numerical results and experimental results.

Theory and method
Biot (1956) considers a volume of the solid-fluid system representation using a cube of unit size. The stress and strain tensor is separated into two parts: solid phase parts and liquid phase parts, the force component acting on the solid parts of each face of the cube is one part. This is denoted as the tensor The scalar S is proportional to the fluid pressure P according to where b is the the fraction of fluid area per unit cross section. And two phase media stress-strain relations are expressed as

Ne
Ae In this form it is easy to express them for arbitrary directions of the coordinate axes by using the invariance of tensor relations. Here, ij d denotes Kronecker symbol function, e ij denotes the solid parts of strain tensor, e, and e denotes the divergence of solid and fluid displacement vector, respectively, A, N denotes elastic parameters of solid skeleton, Q denotes elastic coupling parameters between solid parts and liquid parts, R denotes elastic parameters of liquid parts, scalar ij s and S denote the solid parts and liquid parts of the stress tensor (Yong et al 2006).
The dynamic equilibrium equation of two phase media is also established with simultaneous kinetic equations and dissipative equations (Biot 1956).
where u u u u , , x y z = ( )denotes liquid phase displacement vector, , 11 r 12 r and 22 r denote equivalent density of solid phase, liquid phase and the coupling of two phases, respectively. b denotes dissipative parameters of relative motion between two phases.
With the time derivatives of the variables in equations (2-4) and (2-5), the elastic wave stress equation is established using the particle velocity values to replace the time derivative of displacement, that is v u, , , (2-7), we obtain the elastic wave velocity equation. ).
For 3D two-phase media, the elastic wave first-order velocity and stress equations of Biot theory can be expressed as and v z denote the velocity component of wave propagation in solid phase in the z-, y-and x-directions, respectively. Also, V V , x y and V z denote the velocity component of wave propagation in liquid phase in the z-, y-and x-directions, respectively.
A 3D staggered grid schematic diagram is shown as figure 1, which is divided into many equidistant grids, including main grid and staggered grid. The elastic wave field component and elastic parameters are defined in both grid points, respectively, as table 1 shows (Pei 2005, Pei 2006a, Pei 2006b, and the stress-strain relations are shown in equations (2-4) and (2-4).

In the case of a grid with cuboid elementary cells of length
difference L x of function f can be expressed as follows, with the definition as above: L We can obtain L x + and L x difference scheme with first time derivative: ) Expressions of other wave field components can also be obtained in analogous calculations. With the definition above and equations (2-13) and (2-14), numerical difference operators are obtained that perform the spatial and temporal derivatives in the z-, y-and x-directions in Similarly, we can obtain other difference operators, therefore, each derivative in the z-, y-and x-direction is replaced by a linear combination of several derivative terms. Considering the PML boundary conditions (Jean-Pierre 1994), the staggered-grid high-order difference format can be derived from equations (2-8)-(2-11) with wave field parameters difference operators. Finally, we can obtain the following results for the derivatives in the z-, x-and y-directions x t x x y t y y z t z z where the components of v t+ ( ) in the z-, yand x-directions are represented by v , .
The algorithm in the C programming language needs 97 3D-arrays to represent all these variables in the equation of staggered grid FDM. Take a homogeneous porous model, for example, of size 200 m×200 m× 200 m. Setting the data type to float type, each of the float data takes 8 bytes, so this algorithm requires about 5.78 GB of memory for all the matrices in equation (2-19). It takes 38 557.722 s to perform this model on an ordinary computer, the grid size and time step of which is 10 m and 0.001 s, respectively. Thus, it is suitable for subsequent computer solving.
The way to estimate the effective velocity of both the compressional and the shear waves is also the key point of the new method. The basic idea of using FDM to numerically simulate rock elastic moduli is to model the dynamic wave propagation using this method and study the velocity of the elastic waves based on a 3D digital core in the long wavelength limit (pore size=wavelength) (Saenger et al 2011).
For the purpose of estimating EEP, the digital rock sample A (porous media) is embedded in a homogeneous elastic region, (yellow region) figure 2(a), called the experimental model. Also, another digital rock sample, B (homogeneous media), is embedded in a homogeneous elastic region, (yellow region), figure 2(b), called the reference model. The elastic properties of the grain material in digital rock sample B and the homogeneous elastic region is the same as that of skeletal material in digital rock sample A. The Fontainebleau (see section 3) digital core consists of 0 and 1; 0 represents the pore and 1 represents the framework. In this paper, we take the Fontainebleau as the digital rock sample A, also, we establish digital rock sample B which consists only of 1 based on Fontainebleau by algorithm.
The same plane wave source was applied to generate waves at the top of the two models. With two receivers at the bottom of the two models, it is possible to measure the time delay of the peak amplitude of the plane wave caused by the inhomogeneous structure of the digital rock sample. With the time delay (compared to the reference model) one can estimate the effective velocity of both the compressional and the shear waves (Saenger et al 2011). As a result, we can calculate the effective rock moduli of porous media. The formulas for calculation are as follows.
where S1 denotes the height of the homogeneous elastic region, S2 denotes the height of the 3D digital core sample, v 0 denotes the velocity of the homogeneous elastic region, v eff denotes the equivalent velocity of the digital core sample. A t , 1 t 2 denotes the time the wave propagates through the experimental model and reference model, respectively. where v p and v s denote equivalent compressional and shear wave velocity, respectively, K eff and eff m denote equivalent bulk moduli and shear moduli, respectively; , f r s r and eff r denotes fluid density, framework density and effective density, respectively; j denotes the porosity of the 3D digital core. The equivalent elastic moduli of the 3D digital core can be obtained with equations (2-28), (2-28) (Yin et al 2016).

Model test based on Fontainebleau
This paper establishes 3D digital core models with different porosities utilizing the 3D digital modeling method. For the Fontainebleau sandstone digital core, which has a porosity of 20%, we replace the pore pixels at the edge of the pore with matrix pixels to reduce porosity. Conversely, the porosity can also be increased. In this section, five digital core models with different porosities, 5%, 10%, 15%, 20% and 25%, were constructed (Yin et al 2017). We estimate equivalent elastic moduli with FDM based on these five digital samples. We compare different rock physics models and FDM for the effective elastic moduli, and analyze the equivalence between the digital core results and rock physics theory. Also, we compare staggered-grid FDM numerical results with ordinary-grid FDM numerical results as well as experimental results (Han et al 1986) to test and verify the feasibility and accuracy of the new method.
The size of Fontainebleau is 200×200×200 pixels, and the resolution of Fontainebleau is 5.68 μm/pixel. Fontainebleau is an ideal 3D digital core with framework and pore structure ( figure 3). Fontainebleau sandstone has extensive experimental data regarding elastic moduli, which makes it easy to compare the experimental data and numerical simulation results for various porosities. Fontainebleau's matrix contains only quartz grains and no shale, and its pore structure is only intergranular pores.

Equivalence analysis with theoretical rock physics models
The science of rock physics addresses the relationships between geophysical observations and the underlying physical properties of rocks, such as composition, porosity, and pore fluid content (Mavko et al 1998). For the purpose of analyzing the equivalent conditions between theoretical methods and FDM, we compare the results simulated by FDM and theoretical rock physics models, such as Hashin-Shtrikman bounds theory, SCA theory (Mavko et al 1998) and Nur theory.
Hashin-Shtrikman (H-S) is boundary theory, which represents the narrowest range of elastic moduli values of the isotropic two phase media and can be used to estimate the rock moduli more accurately than Voigt-Reuss bounds (Mavko et al 1998). The bounds of the bulk moduli and shear moduli can be calculated by this theory. Figures 4(a), (b) show the comparison of the Hashin-Shtrikman bounds and the FDM simulation based on the digital cores. In this figure, the bulk moduli and shear moduli of the digital core decrease with the increasing porosity, which is consistent with the trend of the Hashin-Shtrikman upper bounds. The FDM results based on the digital core are distributed within the Hashin-Shtrikman bounds, indicating that the results are reasonable and the FDM can be used to calculate the elastic moduli of the rock.
SCA theory considers the interaction in inclusions and uses the effective medium to replace the background medium. It is an effective medium theory, which is used to calculate the effective moduli of a two-phase medium for high frequency. SCA theory takes the effect of pore shape into consideration, which means that the variation in effective elastic moduli with porosity is related to the pore aspect ratio of the rock. In SCA theory, the porosity is not added to the solid matrix and is instead added to the equivalent medium. The elastic moduli of the equivalent matrix changes constantly until the displacement field generated by the equivalent medium is equal to the displacement field caused by the scattering of the porous rock (Mavko et al 1998). We calculate bulk moduli and shear moduli with different   Figures 4(c), (d) show that, as the porosity of the digital core changes from 0.01 to 0.25, the trend of the FDM simulation is consistent with the trend of SCA theory for pore aspect ratios between 0.2 (the second curve) and 0.40 (the fourth curve). And each of the pore aspect ratios of the bulk moduli in five porosity in figure 4(c) is 0.1 smaller than that of shear moduli in figure 4(d). In particular, when the porosity is from 0.1 to 0.15, and also 0.25, the bulk moduli and shear moduli of the digital core coincide with the third curve (the pore aspect ratio equal to 0.3) and second curve (the pore aspect ratio equal to 0.2).
The theoretical analysis and numerical simulation indicate that, for SCA theory, equivalence is established in cases where the pore aspect ratio is within a certain range. According to the equivalence, the comparison of the SCA theory and DRP for the elastic moduli is a good approach for predicting the pore aspect ratio of the rock.
Nur theory is a dry frame model; it expresses the relationship between the frame moduli and the matrix moduli based on critical porosity (Nur 1992).
Figures 4(e), (f) show the comparison of the Nur theory and the FDM simulation based on the digital cores. In figures 4(e), (f), the trend of FDM simulation based on the digital core is consistent with the trend of the Nur model with the porosity from 0.05 to 0.2, and the former is scattered around the latter. However, it can be seen that the values from both methods are not well matched, which means that these methods are not equivalent.

Comparison of numerical results with experimental results
Han et al (1986) measured compressional velocity V P and shear velocity V s as functions of pressure in 75 sandstone samples, with porosities f ranging from 2 to 30 percent and volume clay content C ranging from 0 to 50 percent, respectively. The velocities are measured at ultrasonic frequencies using the pulse transmission technique (Birch 1960). Han's study investigated the effect of clay content versus porosity on the velocities in sandstones systematically, under laboratory conditions, using a large number of samples covering a wide range of porosities and clay contents. In this section, we calculate the rock moduli according to the experimental value from Han et al (1986). Figures 5, 6 show the comparison between experimental results and the FDM simulation based on the digital cores. In figures 5, 6, the rock moduli (blue rhombus) were estimated using the high-order FDM simulation in a staggered-grid based on the digital cores, which is the new method introduced in this paper. In order to compare the simulation using the new method with the FDM method, we take the rock moduli (gray triangular) which were estimated by the FDM in ordinary-grid (Qin 2017) into consideration. The experimental results (orange squares) are the rock moduli calculated by equations (2)-(2) with the value of velocity and porosity from Han et al (1986).
In figure 5, as porosity increases, we notice that the bulk moduli (staggered-grid FDM) agree well with experiment results. In particular, when the porosity is 5% and 10%, the two are almost identical. When the porosity is greater than 10%, the bulk moduli (staggered-grid FDM) are slightly less than the experimental results. In figure 6, as porosity increases, the shear moduli (staggered-grid FDM) are consistent with experimental results in the porosity of 10%, 15%, 20%, which agrees better than expected for bulk moduli with 5% porosity. In general, the experimental results were consistent with staggered-grid FDM based on DRP when porosity was 10%, 15%, and 20%, which verifies the accuracy and feasibility of this method.
Compared with the rock moduli estimated by the FDM in ordinary-grid, it is easily noticed that the bulk moduli (staggered-grid FDM) are more consistent with experimental results than the bulk moduli (ordinary-grid FDM) for every porosity (figure 5). When the porosity is 5% and 20%, the Figure 5. Comparison between the bulk moduli estimated by FDM and experimental data. Blue rhombuses represent the bulk moduli estimated by the staggered-grid high-order FDM, gray triangles represent the bulk moduli estimated by the ordinary-grid high-order FDM, orange squares represent the bulk moduli estimated by the experimental data. Figure 6. Comparison between the shear moduli estimated by FDM and experimental data. Blue rhombuses represent the shear moduli estimated by the staggered-grid high-order FDM, gray triangles represent the shear moduli estimated by the ordinary-grid high-order FDM, orange squares represent the shear moduli estimated by the experimental data. difference between them is more obvious. In figure 6, when the porosity is 5%, the shear moduli estimated by ordinarygrid FDM is nearer to experimental results than that which was estimated by staggered-grid FDM. However, the opposite conclusion can be drawn when comparing the two methods with porosities of 10%, 15% and 20%, with slight differences. Therefore, rock moduli estimated by staggered-grid FDM is more accurate than those estimated by ordinary-grid FDM, especially in porosities of 10%,15% and 20%.
According to the analysis, the comparison of the experimental results with staggered-grid FDM based on DRP, we have tested and verified the accuracy and feasibility of the staggered-grid FDM, which is the new method. However, the numerical results from the new method still deviate from experimental results in several cases, which prompts us to continue improving the precision of the algorithm. The comparison of ordinary-grid FDM with staggered-grid FDM based on DRP shows that the results of staggered-grid FDM are nearer to the experimental results in a certain range of porosity. Thus, staggered-grid FDM is an effective method to estimate rock moduli based on DRP.

Conclusion
In this paper, we perform the staggered-grid high-order FDM to estimate the effective elastic moduli of two-phase media based on the Fontainebleau digital core, which is easy to implement and more accurate than theoretical rock physics models. This paper examines the equivalence between the numerical results from this new method and various rock physics models, and thoroughly analyzes the conditions of the equivalence relationships. The staggered-grid FDM simulation based on DRP is equivalent to the Hashin-Shtrikman bounds theory, and the DRP model is equivalent to the SCA theory in cases when the pore aspect ratio is within a certain range. According to the equivalence relationships, we are able to predict the pore aspect ratio of the rock. Meanwhile, the numerical results estimated by the new method are consistent with experimental results in a certain range of porosities, and are also more accurate and more feasible than those estimated by ordinary-grid FDM. Thus, the new method of staggeredgrid high order FDM is an effective method to estimate rock moduli based on DRP; it may also be a new method for estimating the other EEP of media numerically, which should be studied further.