1 Introduction

The wave propagation in heterogeneous porous media is one of the long standing topics in geophysics. The fundamental challenges arise mainly from two aspects, i.e., multiphase materials and broad distribution of microstructures. Gassmann^[1] proposed the fluid substitution model, treating the elastic moduli and wave velocities in the uniform porous rocks fully saturated by fluids. Biot^[2-4] proposed a well-known poroelasticity theory to describe the elastic wave propagation in the porous media saturated by viscous fluids. Gassmann's formulas are actually the lower frequency limit of Biot's equations. Because of the complexity of the microstructures of real-world rocks, the pore compliances are not uniform in many cases, leading to unbalanced fluid pressure and local fluid flow among different regions at the mesoscopic scale. A single effective macroscopic fluid flow is assumed in Biot's theory. However, consistent evidences have shown that the attenuations related to the heterogeneous distribution of the fluid saturation have been underestimated by Biot's equations^[5-9]. In order to account for the observed attenuations in seismic data, many efforts have been made to develop other energy dissipation mechanisms, such as the patchy saturation model^{[5, 10-15]} and the squirt flow model^{[8, 16-17]}.

Patchy saturation models usually concern spatial fluid patches, which are solely saturated by immiscible fluids. The high compressibility contrast between fluid patches produces significant attenuation and velocity dispersion. During the last two decades, there are more and more endeavors in understanding the effects of the matrix heterogeneities on velocity dispersion and attenuation^{[5, 18-20]}. Berryman and Wang^[21] and Pride and Berryman^[22] extended Biot's poroelasticity by making a generalization for a double-porosity/dual-permeability model, and analyzed the poroelasticity to seismic wave propagation^[23]. Pride et al.^[5] studied the wave induced flow due to the heterogeneity in the elastic moduli. Müller and Gurevich^[24-26] studied the impact of spatial permeability fluctuations, and provided a dynamic-equivalent-permeability model to predict the P-wave attenuation. Rubino and Holliger^[20], Carcione and Picotti^[27], and Ba et al.^[28-29] studied the effects of pore fluid and solid frame heterogeneities on the dispersion/attenuation curves. Quintal et al.^[30] reported a finite element model for the seismic attenuation/dispersion caused by the compressibility difference. Sun et al.^[31] compared some of the main patchy models by analyzing low-and high-frequency asymptotic behaviors, and used them to experimentally measured the velocity-saturation relations.

The porous media containing patchy saturations and multiphase matrix components can be regarded as a continuum consisting of various composite materials. In the theories of composite materials, the analytical determinations of the bulk properties have been studied for more than one hundred years. A composite material is usually regarded as statistically isotropic when the effective constitutive relation is independent of the coordinate system. One of the primary models is composed of spherical/ellipsoidal particles (Material 1) embedded in matrix (Material 2) with a dilute concentration^[32]. Hashin^[33] and Hill^[34] studied the finite inclusion concentration. Brown^[35] and de Loor^[36] showed that the overall bulk property of the composite material was not completely determined by the volume fractions and properties of the phases. The presence of microstructure heterogeneity makes it quite complicated to derive an expression for the exact effective modulus. Hashin and Shtrikman^[37-38] converted the efforts to the determination of the upper and lower bounds of such quantities. A well-known effective medium theory is the self consistent scheme (SCS), which can estimate the average strain in particle phases^{[32, 39-42]}. Voigt^[43], Reuss^[44], and Hill^[45] studied the effective elastic moduli of the composite materials containing polycrystalline aggregates. The Voigt-Reuss-Hill bounds have been improved by many works^{[37, 46-48]}. The poroelasticity of the bodies containing small scale cracks and pores attracts continuous interests^[49-55]. Because of the complexity of the microstructures of real-world composite materials, the derivation of the effective modulus directly from the precise geometric description of microstructures still faces many practical difficulties. The mechanism of how microstructure governs the wave propagation in a heterogeneous porous medium is still unknown. An exhaustive review of these topics is beyond the scope of this paper. The comprehensive literature survey of these subjects can be found in Refs. ^[56] and ^[57].

This article will discuss a method deriving the elastic moduli of the heterogeneous media consisting of bimaterial matrices and double-concentration inclusions. Unlike many other efforts assuming the uniform inclusion concentration throughout the single-material matrix, this work tries to develop a heterogeneous matrix model with double porosities. The whole porous material is represented by a volume containing two mixed skeleton materials, i.e., the host matrix and the patchy matrix. The patchy matrix is regarded as the embedded patches made of materials (such as kerogen) different from the continuous host matrix (such as mineral). In addition, both the matrices contain particle inclusions at different concentrations. Following Hashin's work^[58], we derive the bulk moduli for the whole system with a multi-level effective modulus method. By setting the elastic constant and density of the particles to be zeros, the moduli of the bimaterial matrix with double inclusion concentrations are calculated. The newly derived modulus formulas are used to predict the kerogen content and total organic carbon with the help of the wave velocity measured in experiments. Pretty good total organic carbon predictions are obtained for 16 black shale samples^[59] and 21 tight-oil rock samples.

2 Composite material containing bimaterial matrix and double inclusion concentrations

We propose a method to estimate the moduli of the composite material containing a bimaterial matrix and double inclusion concentrations (BM-DIC). Consider a body containing a large amount of non-uniformly distributed inclusions (see Fig. 1(a)). The body is divided into representative cells composed of one spherical inclusion embedded in the matrix material. The blank cells and shaded cells are made of different matrix materials (see Fig. 1(a)). The nth cell volume V_n is enclosed by a surface S_n. The volume concentration of the inclusion (VCI) of the nth cell is defined by

where v_n is the volume of the nth inclusion, and a_n and b_n are the radii of the inclusion and the concentric matrix shell, respectively (see Fig. 1(c)). The overall VCI is

where V is the volume of the whole body. The shaded cells have different inclusion concentrations compared with the blank cells in this derivation. Since the inclusions are non-uniformly distributed in the body, the VCI within the arbitrary volume element (much larger than the inclusion size) is not expected to be close enough to c.

The surface S_n can be constructed in any irregular shape. Nevertheless, a spherical surface will be chosen to approach S_n. One typical cell volume is regarded as concentric spheres with the radii a_n and b_n (see Fig. 1(b)).

Let a constant isotropic radial stress

be applied to the body surface, which has linear displacements, the corresponding boundary value problem of which has been solved by Love^[60]. σ⁽⁰⁾ is given by

where σ_ij⁽⁰⁾ are stresses in an elastic homogenous and isotropic body without inclusions. The stresses inside the inclusion are^[60-62]

(1)

(2)

where K and G are the bulk and shear moduli, respectively. The subscripts p and m indicate the moduli of the inclusion and the matrix materials, respectively. Substituting σ⁽⁰⁾ and σ^(p) into the elastic energy perturbation caused by the nth inclusions, we have

(3)

Based on the theorem of the minimum complementary/potential energy, the actual state of the stress minimizes the strain energy U^{(σ)[58, 60]}. Then, we have K^*≥K₁^*, where

(4)

If the isotropic radial displacement is prescribed on the body surface, the energy perturbation δU_n^(ε) cause by the nth inclusion can be obtained in an analogous way as δU_n^(σ). According to the theorem of the minimum complementary/potential energy, the actual displacements minimize the strain energy. Then, we have

i.e.,

(5)

Although K_l^* and K_u^* provide approximated bounds for the bulk modulus, it is to be expected that the relationship between the inclusion modulus and the matrix modulus should be derived so as to guarantee the bounds coincidence in the materials with heterogeneous inclusion concentrations.

3 Methods

Consider a body containing N inclusions, which have two types of concentrations, i.e., low concentration c_l and high concentration c_h. The matrix materials at low and high concentration regions are not necessarily the same. Denote the moduli of the low concentration matrix (the host matrix) by (K_m, G_m), the moduli of the high concentration matrix (the patchy matrix) by (K_m^h, G_m^h), and the moduli of the inclusion materials for the host and patchy matrices by (K_p, G_p) and (K_p^h, G_p^h), respectively. The number of the inclusions with the concentrations c_l and c_h are N_l and N -N_l. In order to simplify the derivation, we choose the same cell radius b_n=b for all the composite cells. Then, we have the volume fraction

We consider the patchy matrix with the high concentration inclusion aggregation as a single effective "inclusion" with the moduli (K_p, G_p) (see Fig. 1(b)). The effective single-inclusion contains N_h particles with the moduli (K_p^h, G_p^h). The particles are surrounded by the matrix material with the moduli (K_m^h, G_m^h). The effective moduli (K_p, G_p) of the patchy matrix are obtained by^[58]

(6)

(7)

where ν_m^h is Poisson's ratio of the patchy solid matrix. The effective host matrix moduli (K_m, G_m) are derived by^[58]

(8)

(9)

where ν_m is Poisson's ratio of the host matrix. It is clear that when the bulk moduli K_m, K_p are the same as (K_m^h, K_p^h) and the concentrations c_l and c_h are equal, the effective moduli K_p, G_p of the heterogeneous region are exactly the same as the background effective moduli (K_m, G_m). This is the uniform concentration condition adopted by Hashin^[58]. When the high concentration c_h region is made of different materials or c_h is quite different from c_l, the overall elastic constant of the composite material can be estimated by the bounds (K₁^*, K_u^*) according to K^*=K₁^*=K_u^*.

(10)

(11)

where α is the fraction of the patchy matrix defined by

For the shear modulus, the up and low boundaries are not overlapped with each other. When the boundaries come close, Hashin^[58] suggested a formula of effective shear moduli as follows:

(12)

where ν_m is Poisson's ratio of the effective host matrix defined by

i.e., the overall Poisson's ratio of the host matrix with inclusions.

4 Examples 4.1 Carbonate rock with double inclusion concentrations

Here, we use the BM-DIC model to calculate the moduli of the carbonate rock containing double inclusion concentrations. The predicted modulus-concentration curves clearly show that the BM-DIC model is consistent with the Hashin model, and non-uniform inclusion concentrations bring entirely different modulus bounds for the composite material.

The inclusion parameters K_p, G_p and K_p^h, G_p^h are for water. Since the fluid cannot bear tensile and shear loads, the shear modulus of water remains zero. The calculated effective bulk moduli of the rocks correspond to the compression case. The host/patchy matrix moduli are from carbonate. The inclusion concentrations in the host and patchy matrices are c and c_h, respectively. The rock parameters are listed in Table 1.

We calculate the cases of the uniform concentration (c_h=c) and the heterogeneous concentration (c_h=2c). In the first case, the BM-DIC model degenerates to a homogeneous porous medium by using the uniform inclusion concentration c=c_h and the homogeneous matrix material with the moduli K_m=K_m^h and G_m=G_m^h. The results in Fig. 2 show that, the effective moduli of the host matrix (K_m, G_m) are exactly the same as those of the patchy matrix (K_p, G_p) (see Fig. 2(a)), and the bounds of the bulk moduli of the whole composite material predicted by the BM-DIC model are the same as the Hashin bounds. Both the BM-DIC model and the Hashin model are between the Voigt-Reuss bounds (see Fig. 2(b)).

In the second case, the inclusion concentration c_h in the patchy matrix is set to be 2c, where c is the inclusion concentration in the host matrix. The material in the host and patchy matrices are the same. The curves in Fig. 3 indicate that, the effective moduli of the host matrix are larger than the patchy region (see Fig. 3(a)). This is caused by the lower inclusion concentration in the host matrix. The existence of more water inclusions weakens the patchy matrix. Figure 3 also shows that the modulus bounds of the BM-DIC model are lower than the Hashin bounds (see Fig. 3(b)). The reason of this lower bounds is that the double concentrations (c_h and c) are used in the BM-DIC model. The higher water concentration c_h=2c causes lower BM-DIC bounds than the Hashin model. In addition, both the BM-DIC bounds and the Hashin bounds are inside the Voigt-Reuss bounds.

4.2 Carbonate/sandstone-clay matrices with double inclusion concentrations

In the second example, the bulk moduli are predicted for a bimaterial matrix containing water inclusions. The host matrix is carbonate, while the patchy matrix is sandstone-clay. The uniform inclusion concentration case (c_h=c) and the non-uniform inclusion concentration case (c_h=2c) are considered, respectively. The primary importance of the predicted moduli-concentration curves is in that the heterogeneity in the matrix material can bring non-negligible deviations from the uniform matrix model. The dispersed patchy matrix regions will weaken or strengthen the bulk moduli of the whole composite material, which is not handled by the Hashin model. The rock parameters are listed in Table 2.

The predicted bulk moduli in Fig. 4 imply that, the effective bulk moduli of the host matrix are larger than the patchy region. This effect is reasonably caused by the weaker material (sandstone-clay) in the patchy matrix, even if the inclusion concentrations are the same, i.e., c_h=c. Figure 4 also shows that the effective shear moduli are nearly the same as the patchy region, and the bounds of the BM-DIC model are lower than the Hashin bounds. This is also the effect of the existence of the softer patchy matrix material.

Note that the Hashin bounds are even higher than the Voigt bound when the inclusion concentration c is less than 0.2. These unexpected cross arises from the fact that we only use the host matrix parameters to calculate the Hashin bounds. This means that the whole porous medium is regarded as a uniform carbonate porous rock saturated by water. While in calculating the Voigt bound, both carbonate and sandstone-clay matrix materials are used. The bounds in such cases are no doubt lower than the Hashin bounds.

In another case, we introduce the double water concentrations in the BM-DIC model. The patchy matrix has a higher water inclusion concentration c_h=2c. Figure 5 shows that the effective bulk modulus of the host matrix is larger than the patchy region. This is caused by the soft material and higher inclusion concentration in the patchy matrix. The effective shear modulus in the host matrix is larger than the patchy region, but nearly the same when c_h is near 0. This is obviously caused by the water concentration in the patchy matrix region. The bounds of the BM-DIC media are lower than the Hashin bounds. The reason of the lower BM-DIC bounds and the unexpected higher Hashin prediction is the same as that in the uniform concentration case.

The results show that the heterogeneity of the matrix materials and porosity in the BM-DIC model can lead to considerable differences in calculating the elastic moduli. Unlike the conventional matrix-inclusion models, the BM-DIC provides extra variables, which can improve the prediction of the elastic moduli when there are composite mixtures in the reservoir rocks such as various mineral contents and kerogen.

In reality, the pores filled with fluids in the reservoir rock are usually under high pressure. For the sake of simplification of the derivations, the effects of the initial fluid pressure on the effective properties of the porous rocks will not be taken into account here.

4.3 Kerogen content estimation in black shale

Here, we use the BM-DIC model to calculate the kerogen content of the Mississippian-Devonian Bakken black shale^[59]. The ultrasonic velocities at different depth follow those in Ref. ^[59]. The velocities are along the normal direction to bedding (θ=0). The inorganic matrix grain moduli are those of clay (K_m=25 GPa, G_m=9 GPa). The kerogen moduli (K_p=4 GPa, G_p=1 GPa) follow Ref. ^[62]. According to the lab data^[59], we give a variation interval (0.001, 0.45) to the patchy (organic material) matrix volume fraction, which is regarded as the kerogen content.

The simulated annealing (SA) method is used to inverse the patchy matrix volume fraction by minimizing the difference between the calculated and the measured bulk moduli. The seven parameters to be optimized are (i) the elastic moduli K_m, G_m, K_m^h, and G_m^h; (ii) the host (mineral) matrix porosity c and the patchy (organic) matrix porosity c_h; and (iii) the patchy (organic) matrix volume fraction α. The initial values for the moduli are those of clay and kerogen. The matrix elastic moduli have a varying range [50 %, 150 %]. The initial porosity and density of the rock samples have been given by Vernik and Nur^[59]. The bulk moduli and density of the particle inclusions are set to be zeros to simulate the pore space.

The inversed kerogen volume fraction α are transferred to the values of the total organic carbon by

where ρ_krg and ρ_b are the kerogen (patchy matrix) and the bulk (host matrix) density, respectively. The constant a=0.75 relates the kerogen weight percent to the total organic carbon T_oc^[59]. The inversed T_oc for the 16 rock samples are shown in Fig. 6(b). The inversed bulk moduli of the samples are compared with the lab values (see Fig. 6(a)).

We find that the calculated bulk moduli match the lab values (calculate from the P/S wave velocities and the bulk density) very well for most samples. The values of the inversed T_oc are consistent with the measured values^[59], although there are misfits for each sample. The optimization errors arise from several aspects. Firstly, the measured velocities are for 1 MHz P wave and 0.7 MHz S wave. Wave dispersion may lead to a decrease in the velocity at a low frequency limit. In addition, the overall measurement error is about 1 percent. Secondly, the kerogen-rich shale samples contain microcracks and thin layers. Strong velocity anisotropy can be observed in different directions^[59]. Here, we use the velocities measured along the bedding-normal symmetry axis. The bulk modulus variations with respect to the bedding-normal directions may lead to uncertainties in predicting the values of T_oc. An anisotropic TM-DIC model is to be developed in a future article.

4.4 T_oc prediction for tight-oil rocks

In this example, we predict the values of T_oc of the tight-oil rock samples which are composed mainly of siltstone, dolomite, limestone, and mud rocks. The rock samples are retrieved from the depths of about 3 100 m to 3 300 m. The values of T_oc and the ultrasonic P-wave/S-wave velocities have been measured on these samples. The pulse transmission technique is used in the velocity measurements with 500 kHz central frequency for the P-wave and S-wave transducers. All the velocities are measured on the dry core samples over a range of confining pressure (5 MPa to 40 MPa). The overall accuracy of the velocity measurements is about 1 %. The porosity is measured at the ambient conditions with a helium porosimeter. Figure 7 shows the measurements of the P-wave and S-wave velocities, densities, and porosities.

We use seven parameters as the inputs to calculate the overall effective bulk modulus by the BM-DIC model. These parameters are the bulk and shear moduli of the host/patchy matrices K_m, G_m, K_m^h, and G_m^h, the host (mineral) matrix porosity c, the patchy (organic material) matrix porosity c_h, and the patchy matrix volume fraction α. The patchy matrix volume fraction α is inversed as the organic material (kerogen) content by the SA method. The moduli of the particle inclusions in the host/patchy matrices are set to be zeros to simulate the pore spaces. The measured porosity are regarded as the initial values of the pore concentrations in the host/patchy matrices for the global optimization problem.

The objective function of the SA is the difference between the predicted bulk modulus K_b^inv and the laboratory bulk modulus K_b^lab calculated from the measurements of the density and the P/S wave velocities. The variation range of α is from 0.001 to 0.45. The initial values for the host matrix moduli are calculated from the velocity and density measurements. The initial values of the patchy matrix moduli are those of kerogen (K_p=4 GPa, G_p=1 GPa)^[62]. All the elastic moduli have a variation range [50 %, 150 %] for the SA algorithm.

The predicted T_oc and experimental measurements are shown in Fig. 8. Since there are seven variables to be optimized and the objective function is based on the complex formulas of the bulk moduli, the predication of T_oc is a complicated global optimization problem with a large amount of local minimums. The objective function surface is so rugged that even the simulated annealing algorithm cannot converge the same optimized values all the time. Variations have been observed in the predicted values of T_oc. Instead of using the solution of a single SA optimization, the average value of 5 SA runs is regarded as the final prediction.

Figure 8(b) shows that the consistence between the predicted T_OC values and the experimental measurements is pretty good for most of the 21 core samples. The difference between the predicted K_b^inv and the measured K_b^lab, i.e., the objective function values, is very close to zero for most of the rock samples. The results clearly demonstrate that the method in this article is capable to predict the bulk moduli of a complex porous medium containing bimaterial matrices and non-uniform inclusion concentrations. Since the newly induced patchy (organic material) matrix parameters reflect the matrix heterogeneity, the organic content is to be obtained in a natural way. The most promising fact is that when the organic matrix volume fraction approaches the actual kerogen content in the rock samples observed in laboratory, the bulk moduli calculated by the BM-DIC model is quite close to the measured values. Nevertheless, the observable deviations from the laboratory measurements still exist in the prediction of several core samples. Some predicted values are even far from the experimental measurements. It is no doubt that delicate improvements are required for the BM-DIC model, including the incorporation of anisotropy, the wave-induced fluid flow effect, etc.

One of the main advantages of the BM-DIC model is that the kerogen content can be straightforwardly represented by the patchy matrix volume fraction. Minimizing the difference between the predicted moduli and the measured bulk moduli, we can easily convert the optimized kerogen content to T_oc. Thus, we now have a method to predict T_oc in rock samples directly from the measurements of the P/S wave velocities, densities, and porosities.

5 Conclusions

The elastic modulus formulas for a porous medium containing double porosity and double matrix materials are derived. The existence of the patchy (organic material) matrix brings a natural way to incorporate the influence of the distributed kerogen in reservoir rocks. An entirely different assumption from other works is that, in the BM-DIC model, the particle inclusion concentration (porosity) is non-uniform throughout the medium. One of the main advantages of the BM-DIC model is that all required parameters have distinct physical meanings and can be easily measured in laboratory. The BM-DIC model makes it possible to characterize the heterogeneity in skeleton materials and porosities simultaneously. The method is used to inverse the kerogen content in black shale and tight-oil rock samples. A promising fact is that the predicted values of T_OC are in pretty good consistence with the experimentally measured data for most of the rock samples, which confirms the importance of the heterogeneity in the matrix material and porosity for predicting the elastic moduli of reservoir rocks.