Validation of Digital Rock Physics Algorithms

Abstract

With a detailed microscopic image of a rock sample, one can determine the corresponding 3-D grain geometry, forming a basis to calculate the elastic properties numerically. The issues which arise in such a calculation include those associated with image resolution, the registration of the digital numerical grid with the digital image, and grain anisotropy. Further, there is a need to validate the numerical calculation via experiment or theory. Because of the geometrical complexity of the rock, the best theoretical test employs the Hashin–Shtrikman result that, for an aggregate of two isotropic components with equal shear moduli, the bulk modulus is uniquely determined, independent of the micro-geometry. Similarly, for an aggregate of two isotropic components with a certain combination of elastic moduli defined herein, the Hashin–Shtrikman formulae give a unique result for the shear modulus, independent of the micro-geometry. For a porous, saturated rock, the solid incompressibility may be calculated via an “unjacketed” test, independent of the micro-geometry. Any numerical algorithm proposed for digital rock physics computation should be validated by successfully confirming these theoretical predictions. Using these tests, we validate a previously published staggered-grid finite difference damped time-stepping algorithm to calculate the static properties of digital rock models.

1. Introduction

It is now technically possible to construct a high-resolution 3-D image of a real rock sample, and to use this to compute the elastic properties of that rock sample, numerically (e.g., in [1]). The present work examines some of the issues which can arise in this program. Some general remarks on common issues are followed by defining procedures for validating any proposed numerical algorithm, and an application to one particular algorithm.

1.1. Common Issues

Any digital image of a rock sample is defined on a discrete imaging grid, and the subsequent numerical calculation is also conducted on a discrete numerical grid. Typically, the imaging grid is Cartesian, and very fine. The numerical grid is also typically Cartesian, and may be coarser, to save computation effort. However, the numerical resolution should be sufficiently fine that the contents of each numerical cell may be treated as uniform.

Furthermore, it is known that very tiny geometric features, like cracks and grain-contacts, can have a disproportionate effect on the overall elastic response, depending upon the micro-geometry. For example, the elastic effect of a crack depends, among other factors, on its ratio of thickness/length [2]. It is clearly important that the imaging grid, and the numerical grid, have sufficient spatial resolution to fully describe such features, in the sense described above [3].

Then, elastic properties are assigned to the grains. Of course, this means that, in addition to grain geometry, the mineralogy for each grain in the sample image must be determined, this may or may not be straightforward, depending on the sample.

Further, most minerals are anisotropic, and this creates another fundamental difficulty. For some minerals, the internal crystalline symmetry influences the external shape; for others not so much. For example, clay minerals often have the external shape of platelets. For such minerals, assigning the appropriate stiffness tensor and orientation to the grains is straight-forward.

For other minerals, the external shape of a grain may not indicate the orientation of its crystalline axes. In such a case, it may not be clear (without X-ray diffraction studies) how the anisotropic axes of any particular grain are oriented. In general, in some places within the rock, stiff axes will be juxtaposed against soft faces, but not in others. Therefore, the rock will be effectively heterogeneous, in the elastic sense, even if it were macro-homogeneous and mono-mineralic. This may be called “orientational heterogeneity”, as opposed to “compositional heterogeneity”; most real rocks have both.

In an energetic depositional environment, the orientations of the grains will be partially or fully randomized, and so it is tempting to assume that the orientation of the anisotropic axes of such minerals is random. It may be that, if the micro-anisotropy is indeed randomly oriented, the anisotropic stress tensor in each grain can be replaced by its isotropic average, but this has not been demonstrated, either theoretically or numerically.

Any digital rock physics calculation should address these issues of resolution and anisotropy, which may affect its accuracy, even when the micro-geometry is realistic.

1.2. Measures of Isotropic Elasticity

In a physical test on a real rock, the stress is applied externally, and the deformation is measured also on the exterior surface. Hence the incompressibility can be measured from the surface average of interior compressional stress/strain for a cubic sample with edge L as:

$$K_{Surf}=-p{\langle \frac{normaldisplacement}{L}\rangle }_{surf}^{-1}$$

(1)

where p is the applied pressure, and ${\langle \rangle }_{surf}$ indicates a surface average. Similarly, the shear modulus can be measured from the surface average of interior shear stress/strain as:

$${\mu }_{Surf}={\tau }_{shear}{\langle \frac{sheardisplacement}{L}\rangle }_{surf}^{-1}$$

(2)

where ${\tau }_{shear}$ is the applied shear stress. These same averages can be calculated numerically from a digital rock image, recognizing the issues discussed above.

However, from that same image, the volume average of interior incompressibility may also be calculated as:

$$K_{Vol}=-\frac{{\langle pressure\rangle }_{vol}}{{\langle dilatation\rangle }_{vol}}$$

(3)

where ${\langle \rangle }_{vol}$ indicates an average over the volume. Similarly, the volume average of interior shear modulus may be calculated as:

$${\mu }_{Vol}\equiv \frac{{\langle shearstress\rangle }_{vol}}{{\langle shearstrain\rangle }_{vol}}$$

(4)

These volume averages are different in principle from the surface averages, and may be numerically different, depending on the accuracy of the approximations discussed above.

2. Methods

2.1. Validating A Digital Algorithm for Non-Porous Rocks

In view of the several approximations mentioned above, and others not mentioned here (since they are dependent on particular techniques), it is important to validate any numerical algorithm against a theoretical solution. For application to certain elastic problems with simple geometry (e.g., a borehole in a uniform stressed medium), theoretical solutions are available. However, for a realistic rock-physics problem, the complexity of geometry precludes any general theoretical solution.

It is possible to construct a rock physics model with an idealized micro-geometry for which a theoretical solution exists (e.g., for ellipsoidal cracks [3]). However, this is not useful for a realistic rock. The overall elasticity of a realistic rock is, in general, not theoretically calculable because the microgeometry is so complicated. Instead, theoretical upper and lower bounds on the bulk elastic moduli have been proposed. Those derived by Hill [4] are not absolute bounds, since they depend upon a micro-geometric assumption (non-correlation of stress and strain, at the micro level), which might not be true.

However, strict bounds on the overall elasticity were derived by Hashin and Shtrikman [5], independent of these issues, and furthermore are closer together than are the Hill “bounds”. For solid composites of two isotropic constituents (with bulk moduli K_i and shear moduli μ_i and volume fractions f_i), the lower (−) and upper (+) H–S (Hashin - Shtrikman) bounds for the bulk modulus K and the shear modulus μ are

$$K_{-}\equiv K_1^{}+\frac{f_2}{\frac{1}{K_2^{}-K_1^{}}+\frac{f_1}{K_1^{}+4G_1^{}/3}}\le K\le K_2^{}+\frac{f_1}{\frac{1}{K_1^{}-K_2^{}}+\frac{f_2}{K_2^{}+4G_2^{}/3}}\equiv K_{+}$$

(5)

$${\mu }_{-}\equiv {\mu }_1+\frac{f_2}{\frac{1}{{\mu }_2-{\mu }_1}+\frac{f_1}{H_1}}\le \mu \le {\mu }_2+\frac{f_1}{\frac{1}{{\mu }_1-{\mu }_2}+\frac{f_2}{H_2}}\equiv {\mu }_{+}$$

(6)

In Equation (6), the combination of moduli

$$H_i=\frac{5{\mu }_i(K_i+4{\mu }_i/3)}{2(K_i+2{\mu }_i)}$$

(7)

for each constituent appears in the bounds for the shear modulus. Here, it is assumed (following [5]) that $K_1\le K_2$ and also that ${\mu }_1\le {\mu }_2$; this correlation is common. If it should happen that $K_1\ge K_2$ while ${\mu }_1\le {\mu }_2$, then the roles of $K_{+}$ and $K_{-}$, defined in Equation (5), are reversed. If it should happen that $K_1\le K_2$ while ${\mu }_1\ge {\mu }_2$, then the roles of ${\mu }_{+}$ and ${\mu }_{-}$, defined in Equation (6), are reversed.

The bounds on the bulk modulus (Equation (5)) may be combined to show the difference:

$$K_{+}^{}-K_{-}^{}=\frac{4f_1{f}_2{\left(\Delta K\right)}^2}{3\left(K_1^{}+4{\mu }_1/3+f_1\Delta K\right)\left(K_2^{}+4{\mu }_2/3-f_2\Delta K\right)}\Delta \mu$$

(8)

where $\Delta K\equiv K_2-K_1$ and $\Delta \mu \equiv {\mu }_2-{\mu }_1$. Hence, the bounds on K coincide exactly in the special case where the two minerals have identical shear modulus (∆μ = 0), in which case the bulk modulus is given exactly by

$$K_{HS}\equiv K_{-}^{*}=K_{+}^{*}=\frac{(4f_1{K}_1^{}\mu +4f_2{K}_2^{}\mu +3K_2^{}{K}_1^{})}{(3f_1{K}_2^{}+3f_2{K}_1^{}+4\mu )}$$

(9)

(where μ is the common shear modulus, of the two constituents and of the aggregate) for any random microgeometry of the two isotropic constituents. (The special case of Equation (9) with μ = 0 yields the Reuss average [4] bulk modulus, the correct result for such a suspension.)

The generalization of Equation (5) for composites of more than two isotropic constituents was given by Hashin and Shtrikman [5]. For such composites, the upper and lower bounds on K coincide, yielding a unique result, if the shear moduli are equal for all constituents.

Of course, few if any physical aggregates have constituents with identical shear moduli. However, although Equation (9) is not a useful result for physical aggregates, it is a very useful result for numerical analyses, since it poses a necessary criterion for the validity of any digital rock physics numerical algorithm: any numerical algorithm which is proposed to calculate the elasticity of a realistic 3D geometry should first be validated by verifying that it satisfies Equation (9), when equal isotropic shear moduli are assigned to the grains; this may be called the “K-test”.

Similarly, the bounds on the shear modulus in Equation (6) may be combined to show the difference:

$${\mu }_{+}^{*}-{\mu }_{-}^{*}=\frac{f_1{f}_2{(\Delta \mu )}^2}{(H_1+f_1\Delta \mu )(H_2-f_2\Delta \mu )}\left[H_2-H_1-\Delta \mu \right]$$

(10)

Hence, the bounds on the shear moduli coincide exactly, in the special case where

$$H_2-{\mu }_2=H_1-{\mu }_1$$

(11)

For composites of more than two isotropic constituents, the upper and lower bounds on μ coincide, yielding a unique result, if the condition of Equation (11) holds for all constituents.

The condition of Equation (11) serves as a calculation of any one of the four constituent moduli, when the other three are specified. For example, it implies that

$$K_2=\frac{4{\mu }_2(H_1-{\mu }_1-2{\mu }_2^{}/3)}{(3{\mu }_2-2H_1+2{\mu }_1))}$$

(12)

For composites of two constituents which confirm Equation (11), the shear modulus is given exactly by

$${\mu }_{HS}={\mu }_{-}={\mu }_{+}=\frac{f_1{\mu }_1{H}_2+f_2{\mu }_2{H}_1}{f_1{H}_1+f_2{H}_2}$$

(13)

independent of the microgeometry of the two isotropic constituents.

Of course, few if any physical aggregates have constituents which confirm Equation (11). However, although Equations (12) and (13) are not a useful result for physical aggregates, they are a very useful result for numerical analyses, since they pose another necessary criterion for the validity of any digital rock physics numerical algorithm: any numerical algorithm which is proposed to calculate the elasticity of a realistic 3D geometry should first be validated by verifying that it satisfies Equation (13) when the constituent moduli are constrained by Equation (11); this may be called the “μ-test”.

2.2. Validating A Digital Algorithm for Porous, Saturated Rocks

Modeling a rock with fluid-filled porosity involves additional considerations. The theory of poro-elasticity (e.g., in [6,7]) requires consideration of the pore fluid properties, the frequency of excitation, and the hydraulic condition (e.g., open or closed) of the rock, all of which are outside of the present scope. However, it also requires consideration of K_S, the average solid modulus of the grains of the rock; this topic can be addressed here.

Love [8] (see especially Sections 121 and 123 (iii, iv)) proved that, for any shape of a homogeneous isotropic solid, its elastic response to a uniform increase of pore pressure (on all of its surfaces, internal and external) is independent of that shape, with each linear dimension decreasing proportionally, so that its shape is preserved. Hence the bulk modulus of the solid shape is the intrinsic bulk modulus of that solid. The proof is valid for any homogeneous isotropic shape, if the pore space is fully connected hydraulically, so that the fluid pressure is uniform.

This theorem may be applied to a representative volume element of a porous rock in a (hydraulically open) “unjacketed compression” experiment, wherein the increase in external pressure on a rock is balanced by an equal increase in internal pore pressure. If the porosity is fully connected on the time scale of the compression, the pore pressure will be uniform throughout, regardless of the complexity of the geometry. In such a test, the solid is compressed on all sides with the same additional pressure, and Love’s theorem applies.

A straightforward extension of Love’s proof to the case of a heterogeneous isotropic solid then concludes that, in an unjacketed test on such a rock with fully connected porosity, the solid modulus is approximately the intrinsic average modulus of the grains, as calculated above.

2.3. A Particular Algorithm

The tests proposed above are the most general conclusions of the present work. Successful passage of these tests should be demonstrated for any digital rock physics project proposed in the future. A failure in any test might be due to a shortcoming in the numerical algorithm applied, or to any of the other issues mentioned in Section 1.1, which may be different in each project. Below, we apply these tests to a particular algorithm, in a context where none of these other issues occur.

A staggered-grid finite-difference damped time-stepping algorithm has been proposed by Lin et al. [9,10] (“LFZ”) for this digital rock physics problem. The LFZ algorithm produces a solution to the static elastic problem by solving the corresponding elasto-dynamic wave equation, with attenuation artificially added to damp out the kinetic energy (while keeping the strain energy) so that the dynamic solution converges to the static solution. The advantage of the algorithm is its computational efficiency: the time-marching scheme does not require the inversion of a large matrix, as in the static approach. The computational effort in such an inversion can be so significant that approximations permitting a 2D inversion of a 3D model have been proposed [11]. Therefore, the computational efficiency of the LFZ algorithm offers a significant advantage.

The method is based upon the first-order hyperbolic system of Virieux [12]: the equation of motion for continua:

$$\rho (\overrightarrow{x})\frac{\partial v_i(\overrightarrow{x},t)}{\partial t}=\frac{\partial {\tau }_{ij}(\overrightarrow{x},t)}{\partial x_j}$$

(14)

and the linear isotropic constitutive relation:

$$\frac{\partial {\tau }_{ij}}{\partial t}=\lambda (\overrightarrow{x}){\delta }_{ij}\frac{\partial v_k}{\partial x_k}+\mu (\overrightarrow{x})\left(\frac{\partial v_i}{\partial x_j}+\frac{\partial v_j}{\partial x_i}\right)$$

(15)

where v_i is the particle velocity vector (a function of space and time), ρ is the density, $\overrightarrow{x}$ is the position vector, and λ = K − 2μ/3 is the Lame parameter. However, LFZ augmented the equation of motion (14) with a velocity-damping term:

$$\rho (\overrightarrow{x})\frac{\partial v_i(\overrightarrow{x},t)}{\partial t}=\frac{\partial {\tau }_{ij}(\overrightarrow{x},t)}{\partial x_j}-d\rho (\overrightarrow{x})v_i(\overrightarrow{x},t)$$

(16)

where the constant artificial damping coefficient d was chosen to efficiently damp out the kinetic energy, for the study of static elastic problems. The physical realism of the damping is immaterial; it is just an artifice to attenuate the kinetic energy. With the inertial terms damped out, the resulting deformation is independent of the material density.

These equations incorporate local heterogeneity of the elastic moduli and density explicitly (hence no need to consider local boundary conditions), although the model analyzed herein (described below) has piecewise uniform elasticity within each numerical cell. The dynamic algorithm (Equations (15) and (16)) is implemented on a staggered Cartesian numerical grid, similar to that proposed by Virieux [12], shown in Figure 1.

The derivatives are discretized by centered finite differences, with the various quantities evaluated at different positions on the staggered grid, as shown in the Figure 1. The stress components are evaluated at t time units; the velocity components are evaluated at (t + 1/2) time units. As the computation approaches its asymptotic (static) limit, these time differences become irrelevant. This pattern of staggering arises naturally from the centering of the finite derivatives.

Where a grid node separates domains of different elasticity, the elastic moduli and density at that node are assigned appropriate average values, depending on the properties of its eight nearest neighbors in three dimensions. This procedure partially smooths the spatial variation of physical properties, with an accuracy which depends upon the resolution of the grid.

3. Results

3.1. Application of the LFZ Algorithm to Non-Porous Rocks

The Hashin–Shtrikman tests are here applied to the LFZ algorithm. In order to avoid problems caused by mismatch between the numerical grid and the model grid, the calculation was performed for the model shown in Figure 2, with cubic grains. Although this is not a realistic model, with this geometry, any inaccuracies in the calculation can be attributed to other features of the algorithm. To ensure uniform application of external stress, the model shown was surrounded by a uniform jacket (four cells thick) with the properties of water (for applied pressure) or of an average solid (for applied shear stress).

In the model, the cells are randomly assigned to one or another of two isotropic constituents. For the H–S “K-test” in Section 2.1, the elastic parameters are set as K₁ = 13.564 GPa, K₂ = 8.564 GPa, μ₁ = μ₂ = 4.586 GPa, f₁ = f₂ = 0.5. For this model, the theoretical H–S value from Equation (9) is K_HS = 10.700 GPa. The surface-average incompressibility from Equation (1) is K_Surf =10.689 GPa (an error of −0.1%); the volume-average incompressibility from Equation (3) is K_Vol = 10.685 GPa (an error of −0.1%).

The calculation was repeated, with the same model, but with a numerical grid much finer (200 × 200 × 200 cells), so that each cubical grain was represented by 64 cells, instead of by 1. The surface-average incompressibility from Equation (1) is K_Surf = 10.691 GPa (an error of −0.1%). The volume-average incompressibility from Equation (2) is K_Vol = 10.690 GPa (an error of −0.1%). In this instance (with the numerical grid aligned with the model), the increased resolution made little difference in accuracy.

Then, with the same model geometry, the cells are assigned, for the H–S “μ-test”, the elastic parameters are K₁ = 8.564 GPa, μ₁ = 3.236, μ₂ = 3.886 GPa, K₂ = 4.012 GPa (from Equation (12)), f₁ = f₂ = 0.5. For this model, the theoretical H–S value from Equation (13) is μ_HS = 3.546 GPa. The surface-average shear modulus from Equation (2) is μ_Surf = 3.544 GPa (precise to three significant figures). The volume-average shear modulus from Equation (4) is μ_Vol = 3.546 GPa (precise to four significant figures). The incompressibility is not uniquely determined with these parameters.

3.2. Application of the LFZ Algorithm to Porous, Saturated Rocks

To numerically test the procedures of Section 2.2, the algorithm of LFZ was applied to the model of Figure 3. The previous model (Figure 2) is penetrated by straight interconnected channels of water (K_water = 2.25 GPa) as depicted; the porosity is about 35%. The figure shows the inhomogeneous case, with solid matrix composed 50/50 of the two minerals specified in the “K-test” above. Again, this model is not realistic, however, with this geometry, any inaccuracies in the calculation can be attributed to other features of the algorithm. Since a shell of water surrounds this model completely, and pressure is applied to the exterior of this fluid shell, the solid grains are exposed on all sides to the same pressure.

For a model like that of Figure 3, but with homogeneous solid (with solid incompressibility K₂ = K₁ = 13.564 GPa), the surface-average incompressibility (simulating an unjacketed experiment) from Equation (1) is K_Surf = 13.543 GPa (an error of −0.2%). The calculated solid-volume-average incompressibility from Equation (2) is K_Vol = 13.564 GPa (precise to five significant figures).

For the inhomogeneous model of Figure 3, the surface-average incompressibility (simulating an unjacketed experiment from Equation (1)) of the solid portion only is K_Surf = 10.642 GPa, (an error of −0.5%). The volume-average incompressibility from Equation (2) of the solid portion only is K_Vol = 10.650 GPa, (an error of −0.5%). These differences are a measure of the accuracy of the extension of Love’s theorem to the case of inhomogeneous solid, as well as the accuracy of the LFZ algorithm.

There is no equivalent process for determining the average shear modulus μ_S of the grains of the porous rock. If the H–S “μ-test” were applied to a composite with one constituent fluid (μ₁ = 0), the required K₂ from Equation (13) would be negative (hence un-physical), so the H–S “μ-test” is not applicable.

4. Discussion

We briefly discussed here some of the important issues related to digital rock physics computations. These include issues related to resolution (of the imaging grid and the numerical grid), to registration of the two grids, and to macro- and micro-anisotropy. A primary issue is the need to validate the numerical algorithm with a theoretical solution. For the complicated geometry of a digital rock, suitable quantitative validation procedures for any numerical computation of isotropic solid rock elasticity were presented here.

H–S theory [5] provides rigorous bounds for the bulk and shear moduli of an isotropic aggregate of isotropic components. It is well known that for two-constituent aggregates, wherein both constituents have the same shear moduli, the bounds on the bulk modulus of the aggregate coincide, so that its bulk modulus is uniquely determined, regardless of the internal geometry. Although this case is vanishingly rare in real rock, this “K-test” provides a useful validation procedure for any numerical calculation of digital solid rock physics.

A second consequence of the H–S bounds [5] is that if the moduli of the two constituents are constrained by the present Equation (11) above, the bounds on the shear modulus of the aggregate coincide, so that its shear modulus is uniquely determined, regardless of the internal geometry. Although this case is vanishingly rare in real rock, this “μ-test” provides a second useful validation procedure for any numerical calculation of digital solid rock physics.

For a porous rock with homogeneous solid, the bulk modulus may be measured (and calculated) under the condition of equal external pressure and pore pressure. A theorem by Love [8] proves that, in this (“unjacketed”) condition, the bulk modulus of the aggregate is that of the solid itself. By extension, for a rock with inhomogeneous solid phase, the bulk modulus is approximately the average modulus of the solid itself. This provides a useful validation for any numerical calculation of digital rock physics for porous rocks.

These validation tests were applied to the LFZ algorithm [9,10]. The algorithm calculates the static modulus as the asymptotic limit of the corresponding elasto-dynamic equations, with artificial damping, thus avoiding the necessity to invert a large numerical matrix. The algorithm smooths the discontinuities in material properties at each grain boundary; the consequent errors are smaller with higher resolution of the numerical grid.

For the non-porous inhomogeneous solid, the algorithm reproduces the H–S results within a fraction of a percent. For the porous rock, with either homogeneous or inhomogeneous solid, the algorithm reproduces the Love result for the solid, within a fraction of a percent. We conclude that the LFZ algorithm is well-validated for this model, by all these tests. These inaccuracies would presumably be larger for a realistic model, with grains not conforming exactly to the numerical grid, but such errors would be different for each model and each algorithm, so that no generalizations are possible for such errors.

More generally, we recommend that these tests should be applied to any numerical algorithm which may be proposed to calculate the elasticity of any rock model. Of course, these tests, although necessary, are insufficient for a complete validation, as they do not consider all potential sources of inaccuracy. Nonetheless, only after an algorithm has passed these necessary validation tests should it be applied to a case with realistic grain parameters.