An Equation for the Bulk Modulus of Composites Derived From the Effective Medium Theory

Nwonodi, Roland Ifeanyi; Dosunmu, Adewale; Okoro, Emmanuel Emeka

doi:10.1115/1.4055628

Abstract

Bulk modulus has wide applications in well engineering, seismic exploration, waste reinjection, and predicting pore pressure in carbonate reservoirs. However, there is no easy way to obtain accurate values for the effective bulk modulus of rocks. Practically, researchers use rigorous, costly, and time-consuming experiments on core samples. But, stress release and changing rock’s environment have affected the accuracy of results. Also, it is impossible to get accurate values of the effective bulk modulus from theory without accounting for the deformation of microcracks in the rock. Existing models do not consider the presence of microcracks because of the inability to define the positions of cracks relative to one another. Thus, earlier studies introduced approximations to define the upper and lower bounds of values. This study aims to overcome this limitation by accounting for the fluids in the microcracks, apart from those in stiff pores. From the product of the surface area and thickness of the fluid in the microcracks, the authors generated proportionality between the volume of fluid and that of the grain and obtained expression for the crack porosity. Then analytical and numerical techniques were applied to obtain models for the effective bulk modulus. The results show that the presence and magnitude of inclusions reduce the effective bulk modulus significantly. This was validated by a finite element analysis (FEA) using the FEATool run in matlab. In addition, higher volume of fluids in the microcracks makes the rate of change of the bulk modulus with the porosity to be higher.

Introduction

Rocks are composites comprising a solid matrix and fluid-filled pore spaces [1–8] that combine to produce effective elastic properties. One such elastic property is the bulk modulus. This bulk modulus is a parameter that is incorporated into geomechanics, particularly in planning and managing well engineering projects. It is also important when carrying out geophysical explorations with seismic data [9], and it defines how the volume of porous rocks changes with external pressure, e.g., during the pressurization of an oil wellbore. This parameter is useful during drilling cuttings and wastewater reinjection [1] for the safe disposal of hazardous materials derived from drilling activities. The study can extend its application to the prediction of the bulk modulus of glassy/crystalline materials that possess porosity. Furthermore, it can find application in predicting pore pressures in carbonate rocks because of the relationship between bulk compressibility (inverse of the bulk modulus) and effective pressure. The bulk modulus can be determined statically using experiments on core samples or dynamically using data from wireline logs plus theoretical equations [10–12]. Engineers and scientists mostly use and rely upon the static bulk modulus since it represents exact measurements. However, the experimental procedure is very costly to carry out (time-consuming and rigorous for some rocks, particularly shale or shaly rocks) and may yield inaccurate results if the experimenter does not ensure proper sample preservation. Changing the sample from its environment of stress (or stress release) can contribute to poor results from the experiments. The values of the individual components of the rock can differ significantly from that of the macroscopic (or bulk) frame depending on critical parameters related to geological factors in-situ.

There is currently the need to establish suitable ways of estimating the bulk modulus, especially when data for its components are available. For example, the conventional layered-rock approach that attempts to generate predictive models does not accurately simulate the static process. During experimental runs, there will be closure and/or the initiation of microcracks, which the conventional models generally do not consider. If microcracks are present in-situ, there will be fluids inside them plus those in the pores of the rock. The conventional models do not account for the extra fluids in microcracks. And this study provides a way to overcome this limitation.

The basis for defining the average property of a mixture can be in terms of conductivity, molar fraction, mass fraction, weight fraction, or volumetric fraction. And the weight fraction of a species may not be equal in magnitude to the volumetric fraction since they do not represent the same quantity. These bases are applicable in the theory of the effective medium, which describes analytically the bulk behavior of composites. The effective medium approximation uses theories to estimate the physical properties of heterogeneous materials from their components. Belyaev and Tyurnev [13] used the theory to predict the electrodynamic properties of a dielectric medium composed of metallic nanoparticles. David and Zimmerman [14] applied the effective medium theory to analyze the properties of porous rocks. Finally, Wang and Pan [15] used it to predict the physical properties of complex multiphase materials. Thus, one can relate the volumetric contribution of the components to the bulk modulus. However, the conventional methods assume that the volumetric fraction of fluids is equivalent to the porosity of the components without considering the fluids inside microcracks. It is interesting to uncover accurate ways to estimate the bulk modulus of reservoir rocks using porosity and the elastic properties of the constituents. Apart from the fluids in the pores, some fluids exist inside microcracks. This is of much interest because one cannot estimate correctly the elastic moduli of rocks without knowledge of the geometric distribution of the constituents, including those of microcracks [4,16]. To this end, an upper bound and a lower bound were established for the bulk modulus whereby subsequent models are expected to yield results between the limits. It can invoke some level of curiosity when one considers the accuracy of the upper and lower bounds themselves, especially when such bounds do not incorporate the contribution of microcracks apart from assuming a linear combination of the layers. Accounting for the fluids in the microcracks can help to overcome the challenge of not defining the deformation of the microcracks in the stress-strain approach. Thus, in this study, the fluids in the microcracks contribute to the total fluids in the rock. This is to improve the results of the effective medium approach for predicting bulk modulus. The study considers how the bulk modulus changes with porosity, and strives to improve the gaps in models used to predict the effective bulk modulus.

The Effective Medium for Layered Rocks

The sedimentary rock is stratified into layers and comprises various combinations of the pore space and solid matrix extending over any direction. Allowing a linear combination of two layers, one can represent their volumes by that of a rectangular prism, which is the following expression:

V = L B h

(1)

where V is the volume of the object, h is the thickness of the object, L is the length of the object, and B is the width of the object. The thickness of the material may be difficult to determine at the microscopic level, but not at the macroscopic level. Over the same length and width, one can use Eq. to represent the volumes of each layer as follows:

V_{1} = h_{1} L B

(2)

V_{2} = h_{2} L B

(3)

where V₁ is the volume of layer 1, V₂ is the volume of layer 2, h₁ is the thickness of layer 1, and h₂ is the thickness of layer 2. The total volume of the layers comprises the sum of the contributions from each layer:

V_{T} = L B (h_{1} + h_{2})

(4)

where V_T is the total volume of all the layers in the composite. The percentage by volume of each layer is the volume of the layer divided by the total volume, thus:

% V_{1} = \frac{h_{1}}{h_{1} + h_{2}}

(5)

% V_{2} = \frac{h_{2}}{h_{1} + h_{2}}

(6)

where

% V_{1}

is the percentage by volume of layer 1 and

% V_{2}

is the percentage by volume of layer 2. The sum of the percentage by volume of each layer equals unity. The volume ratio is the thickness ratio when the unit is linear and has the same length and width. Similarly, if the material has the same length, but different thicknesses and widths, the volume ratio reduces to an area ratio. The volume ratio in an arbitrary rock may not necessarily be equal to the area ratio or length ratio, even though the sum of each ratio equals unity. Thus, the effective bulk modulus of the two-layered-rock can be calculated with the following expression:

K_{m} = % V_{1} K_{1} + % V_{2} K_{2}

(7)

where K_m is the effective bulk modulus, K₁ is the bulk modulus of component 1, and K₂ is the bulk modulus of component 2.

Some Conventional Methods Developed for the Bulk Modulus

Several theories have been proposed for obtaining the bulk modulus based on the fractional contribution of the constituents and their deformations. There is the Voigt model which yields the main upper bound of values established in the oil and gas industries [16]. This model is one of the commonly and easily applied models for obtaining the bulk modulus. The principle behind its formulation is that the individual components are subjected to the same strain while deforming. That is, the components are linearly arranged. Fjær et al. [4] claimed that the model is the most appropriate for rocks at low values of porosities but the simplest non-interacting approximations are well suited at low porosities. However, its major limitation is that it ignores the deformation of microcracks in the rock, as is the case with the other conventional models. These microcracks are open surfaces in the rock that contains some fluids. The mathematical expression of Voigt’s model is the following:

K_{V} = φ K_{f} + (1 - φ) K_{s}

(8)

where K_v is Voigt’s bulk modulus, K_f is the bulk modulus of the fluid component, K_s is the bulk modulus of the solid component, and φ is porosity. When one compares Voigt’s model to a two-component form of Eq. , comprised of pores and solid matrix, one observes that the author equated the volume fraction of the fluid to the porosity. Hill [17], based on the principle of energy conservation showed that the maximum possible value of the bulk modulus of a composite is the Voigt bound.

The Reuss model is another commonly applied one in the industry, and it yields results for lower bounds of the bulk modulus. It is founded on the same principle as Voigt’s model but assumes that the solid undergoes similar stresses to the liquid during deformation. When a rock is being carried by a liquid, there may be hydrostatic uniformity between the constituents, which makes the rock tend to behave (or flow) as the liquid. Thus, the model seems to act towards the strength of the liquid in the mixture, yielding lower bounds of values. The mathematical expression for this model is as follows [4]:

\frac{1}{K_{R}} = \frac{φ}{K_{f}} + \frac{1 - φ}{K_{s}}

(9)

where K_R is the Reuss’ bulk modulus. The same principle of energy conservation implies that the Reuss model is the minimum possible bulk modulus for rocks. It seems one can obtain a form of this equation from the bulk compressibility of the rock and its relationship with the bulk modulus. This bulk compressibility considers the deformation of the liquid and solids combined, and most liquids can be incompressible. The bulk compressibility of rocks relates inversely to the bulk modulus as follows [18–20]:

C_{b} = \frac{1}{K_{b}}

(10)

where C_b is the bulk compressibility and K_b is the bulk modulus of the rock. Allowing the following simple expression for the bulk compressibility as a function of the compressibility of the matrix and pore space:

C_{b} = φ C_{p} + (1 - φ) C_{s}

(11)

where C_p is the pore compressibility and C_s is the matrix compressibility. Then the inverse of the compressibility of the components can be substituted into Eq. to finally arrive at Eq. , assuming the fluid occupies the pores only. From Eq. , the bulk compressibility will tend to a uniform value if the components have approximately equal compressibility.

The arithmetic average of Voigt and Reuss’ models has been introduced to get results that fall within the lower and upper bounds. This is the Voigt–Reuss–Hills (VRH) model [21], which produces results between the bounds depending on the relative magnitude of the prevailing component in the rock. The equation is the following expression:

K_{VRH} = \frac{[φ K_{f} + (1 - φ) K_{s}] [φ K_{s} + (1 - φ) K_{f}] + K_{f} K_{s}}{2 [φ K_{s} + (1 - φ) K_{f}]}

(12)

where K_VRH is the Voigt–Reuss–Hill bulk modulus. The VRH model will yield relatively accurate results when the components have similar bulk modulus but the model does not work well when there is a high elastic contrast between the components.

Other approaches have been introduced to predict the bulk modulus. Some of these models possess theoretical background but may yield inaccurate/unphysical values, as they do not define the rock completely. Other models introduced more variables that make analysis hard or complex but still produce results between the bounds. Fjær et al. [4] suggested the following expression based on the concept of critical porosity:

K_{f r} = K_{s} (1 - \frac{φ}{φ_{c}})

(13)

where φ_c is the critical porosity and K_fr is the bulk modulus of the frame. The authors assumed that for a set of sand grains, there is a maximum porosity beyond which the solids will not be in contact anymore. They called this maximum porosity the critical porosity. The basis for this formulation is that the bulk modulus of the frame equals that of the solid when porosity is zero but disappears when the porosity approaches the critical porosity. The values of the critical porosity were taken as being equal to those in the standard packing of spheres for uniform or random samples. The equation cannot define the bulk modulus of the liquid component when porosity equals unity as it is not explicitly based on the theory of poroelasticity. It yields large distortions at low values of porosities. The conventional Voigt and Reuss models satisfy this physical interpretation when porosity is zero or one. Thus, the concept of critical porosity requires slight modification or a redefinition of what critical porosity stands for.

Hashin and Shtrikman [22] attempted to reduce the width of the conventional Voigt and Reuss models. They introduced an equation that is more complicated to solve since it requires more elastic parameters than the three basic ones considered so far. The basis for the formulation is the minimization of the potential elastic energy of the rock assuming spherical symmetry. Their equation for the upper and lower bounds are the following:

K^{HS +} = K_{s} + \frac{φ}{\frac{1}{K_{f} - K_{s}} + \frac{1 - φ}{K_{s} + 1.333 G_{s}}}

(14)

K^{HS -} = K_{f} + \frac{1 - φ}{(1 / (K_{s} - K_{f})) + (φ / K_{f})}

(15)

where G_s is the shear modulus of the solid component, HS stands for Hashin and Shtrikman, K^HS+ is the upper bound, and K^HS− is the lower bound of the bulk modulus, which approximates Reuss’ model. The shear modulus requires the relationship between elastic moduli and the use of Poisson’s ratio. Based on the relationship between Young’s modulus and bulk modulus, and between Young’s modulus and shear modulus, one can approximate the shear modulus to 0.7816 multiplied by the bulk modulus of the solid. This is for an average value of Poisson’s ratio of about 0.19. If one ignores the shear modulus in the upper bound, one can obtain Reuss’ model.

Methodology

The approach used to achieve the goals of this study is combined empirical and analytical methods. The authors modified the volumetric contribution of the fluids and incorporated it into the theory of the effective medium. The effective medium approximation means that the sum of the product of the given parameter and its volumetric contribution is the effective parameter:

P = \sum_{i = 1}^{n} P_{i} V_{i}

(16)

where P is the effective parameter, P_i is the value of the parameter for component i, V_i is the volumetric contribution of the components with the parametric value, and n is the number of components considered. An arbitrary volume of rock with definite porosity and bulk moduli of the constituents were used in the study. To predict the bulk modulus of a composite comprising solids and pores, a two-component system was chosen. The rock possesses some fluids in the microcracks, which are a part of the system’s fluid. Thus, the effective volume of fluids in the rock comprised the volume of liquids in the pores and those in the microcracks:

V_{f} = V_{p} + V_{c}

(17)

where V_c is the total volume of fluid in the rock, V_c is the volume of fluids in the microcracks of the solid matrix, and V_p is the volume of fluids in the pores. It is the belief that this yields a more representative value for the volume fraction of the fluid component. To develop a well-defined expression for the fluids in the microcracks, the authors used the product of the surface area open for storage and the thickness of the fluid in the opening. The volume of fluids in the microcracks of spherical solids was used to develop proportionality between the volume of fluids and the volume of solids. From the set of equations generated, a limiting condition was used to analytically obtain the parametric constant in terms of porosity. The authors also used numerical analysis to develop an expression for the bulk modulus. However, the elegance of analytical methods favored the study. From the analytical approach, the authors were able to generalize to develop a model that can predict the effective bulk modulus of a unit of rock comprised of different lithology.

To pursue the numerical approach, the average bulk modulus of the conventional models was used to find a root for the parameter following the Newton–Raphson algorithm [23] and the LP Simple Engine of the Excel Solver. The conventional models selected for the averaging include the Voigt, Reuss, V–R–H, and HS models. While the initial guess was obtained graphically. In addition, the finite element analysis (FEA) was used to validate the results using properties from a selected wellbore. To do this, an elemental 3D volume representative of the size of the rock was selected and gridded using the FEATool. This tool runs on matlab. Then a model for predicting the bulk modulus of the grids was developed and a computer experiment was run to obtain the strains and stresses from the in-built linear elasticity program. The basis for the developed model is the strain energy per volume needed to compress the grids frame by the application of external load. A force of 1000 N was applied to the faces of the element, and the force divided by area gave the magnitude of the external load. This load served as the boundary condition on the faces of the element.

Data from the Tertiary basin of the Agbada formation in the Niger Delta were extracted to aid the analysis. To get the porosity at the interval, the authors applied the Raymer–Hunt–Gardner equation using the data from the transit and mineral [1]. While the data from gamma-ray aided in obtaining the Poisson’s ratio for the Tertiary formation.

Model Development for the Bulk Modulus of the Composite Rocks

The reservoir comprises a solid matrix, fluid-filled and connected pore spaces, and fluid-filled microcracks existing as a unit. Using the volumetric contributions of pores and solid matrix as the basis of a two-component system, Eq. gives the following expression for the effective bulk modulus:

K = % V_{f} K_{f} + % V_{s} K_{s}

(18)

where

% V_{f}

is the percentage by volume of the fluid component, K_f is the bulk modulus of the fluid component,

% V_{s}

is the percentage by volume of the solid components, and K_s is the bulk modulus of the solid components. Here, the volumetric factor of the fluid need not be solely dependent on the volume of the pores as is the tradition but is determined by including the fluids in the microcracks. The total volume of fluids in the rock equals those in the pores and those in the microcracks as given in Eq. . Wang et al. [9] implied that the total porosity comprised components of the crack porosity and the stiff porosity while experimenting on rock cores. If all the fluids are perfectly contained in the pores, then the volume of fluids becomes equivalent to the volume of fluids in the pores. The solid grains may be irregular, but spherical shapes are applicable. In addition, the grains will most likely be randomly packed to achieve the best density. Nevertheless, the volume of the thin film of fluid existing in microcracks can be given as the product of the surface area covered by the fluids and the thickness of the fluids. Using spherical grains and extracting the radius of the sphere from its volume, yields the following expression:

V_{film} = 4.8 h_{film} V_{g}^{(2 / 3)}

(19)

where V_film is the volume of the film of fluids in the crack, which may be the product of the volumetric rate into the crack and time, h_film is the thickness of the film of fluids in the crack, and V_g is the volume of spherical grains. Thus, the volume of the crack fluids is proportional to the two-third power of the volume of solids, and the following expression holds for the crack volume:

V_{c} = K_{c} V_{g}^{2 / 3}

(20)

where K_c is the proportionality constant with units in volume raised to one-third. One can intuitively relate this parameter to porosity and solid matrix, or use numerical analysis to obtain a closed value.

To follow a numerical scheme, combine Eq. with Eq. and substitute into the expression for the bulk modulus. The following expression emerges after defining the porosity from micromechanics:

K = (φ + a (1 - φ)^{2 / 3}) K_{f} + (1 - (φ + a (1 - φ)^{2 / 3})) K_{s}

(21)

The dimensionless parameter a = V^−(1/3)K_c relates the volume to the proportionality constant of Eq. (20). With the use of the root-finding tool, one can use a set of local or regional data to find the value of the parameter that makes the function F(a) equal to zero. Experts use numerical methods (NM) to easily find the values of function/parameters that may be very complex, rigorous, or hard to solve. However, one cannot compare NM to the elegance of analytical solutions, especially when such solutions are obtainable. Numerical analysis should not be used to make generalizations in the petroleum industry where geologic settings can vary significantly from one region to another. The exception is when such generalization is localized.

If there is no grain in the system, the porosity equals one and there will be no microcracks to store the fluids. Similarly, if there is no fluid in the system, porosity equals zero and there will be no fluid in the microcracks. Thus, a limiting condition is the presence of a fluid-phase and solid-phase in the system, and the following equation satisfies the proportional constant of Eq. :

K_{c} = φ (1 - φ) V_{p}^{1 / 3}

(22)

Equation turns into the following expression for the volume of fluids in the cracks:

V_{c} = V_{p}^{1 / 3} V_{g}^{2 / 3} (φ (1 - φ))

(23)

From Eq. (23), the expression for the total or modified volume of fluids becomes the following:

V_{f} = V_{p} + V_{p}^{1 / 3} V_{g}^{2 / 3} (φ (1 - φ))

(24)

Then the expression for the volume ratio of the fluid can be given as follows:

\frac{V_{f}}{V} = \frac{V_{p}}{V} + {(\frac{V_{p}}{V})}^{1 / 3} {(\frac{V_{g}}{V})}^{2 / 3} (φ (1 - φ))

(25)

From micromechanics, the ratio of pore volume to bulk volume is equal to porosity while that of grain volume to bulk volume equals −1 porosity. Thus, the volume ratio of the fluid reduces to the following:

\frac{V_{f}}{V} = φ + (1 - φ)^{5 / 3} φ^{4 / 3}

(26)

The sum of the volumes of the fluids and solids is the total volume of the unit; thus, the volume ratio of the solid is the following:

\frac{V_{s}}{V} = 1 - \frac{V_{f}}{V}

(27)

Therefore, the effective bulk modulus for the composite becomes the following expression:

K = K_{f} (φ + (1 - φ)^{5 / 3} φ^{4 / 3}) + K_{s} (1 - (φ + (1 - φ)^{5 / 3} φ^{4 / 3})

(28)

Following the same line of reasoning for a single unit with two components of solids and fluids, the effective bulk modulus of a n-unit of rocks can be logically expressed as follows:

K = \sum_{i = 1}^{n} K_{f i} (φ_{i} + {(1 - φ_{i})}^{5 / 3} φ_{i}^{4 / 3}) + K_{s i} (1 - φ_{i} - {(1 - φ_{i})}^{5 / 3} φ_{i}^{4 / 3})

(29)

where K_fi is the bulk modulus of the fluid in the ith unit of the rock, K_si is the bulk modulus of the solid in the unit, φ_i is the unit porosity, and n is the number of rock units in the analysis. For example, to analyze the effective bulk modulus of a three-unit rock comprising sandstone, limestone, and shale, the respective bulk moduli of the fluid and solid components of each unit and their porosity needs incorporation into Eq. .

Continuing the development of Eq. requires the use of micromechanics. Following the definition of porosity as the volume of voids per bulk volume, one can define the crack porosity in a similar form, such that the following obtains:

V_{c} = φ_{c} V

(30)

Using the micromechanical definition of porosity and matrix content, the following emerges:

V_{c} = (1 - φ)^{5 / 3} φ^{4 / 3} V

(31)

Thus, comparing Eqs. and indicates that the crack porosity is approximated as the following expression:

φ_{c} = (1 - φ)^{5 / 3} φ^{4 / 3}

(32)

NB: the concept of crack porosity finds application in the enhancement of the permeability of reservoir rocks and geomechanics. This crack porosity (or fracture porosity) is the result of the tectonic fracturing of reservoir rocks, and it is a type of secondary porosity quantifiable using acoustic or resistivity image logs [24,25]. Unfortunately, the image log is not available for every wellbore, and it is costly and hard to estimate fracture porosity. For example, its estimation requires the integration of the frequencies from the image logs and crack aperture into production and stimulation models [25]. Fortunately, the secondary porosity is taken as the difference between the density porosity and the sonic porosity for clastic or non-clastic rocks [24,26]:

φ_{Sec} = φ_{D} - φ_{S}

(33)

where φ_Sec is the secondary porosity, φ_D is the density porosity, and φ_S is the sonic porosity. It is possible to use this difference to validate the model for the crack porosity developed in this study. The conventional models for obtaining formation porosity are the following [1,24]:

φ_{D} = \frac{ρ_{m} - ρ}{ρ_{m} - ρ_{f}}

(34)

φ_{S} = \frac{t - t_{m}}{t_{f} - ρ_{m}} \times \frac{1}{C}

(35)

where ρ_m is the matrix density, ρ_b is the bulk density from the wireline log, ρ_f is the fluid density, t is the transit time observed, t_m is the matric transit time, t_f is the fluid transit time, and C is the compaction factor for correction, which varies between 1 and 1.3. Apart from theoretical models, laboratory analyses are available for the determination of porosity [24,27].

Analysis and Discussions

This section evaluates the effective bulk modulus of a section of the Nembe Well-X in the Tertiary basin of the Niger Delta. This wellbore is likely a future candidate for drill cuttings reinjection, and Fig. 1 shows the data of the wireline log between 5500 ft (1676.95 m) and 6300 ft (1920.87 m). Based on the transit time in the wellbore with a value of 73.4 µs/ft (2640.85 µs/m) and sandstone at 6070 ft (1850.74 m), the porosity is 0.163 using the Raymer–Hunt–Gardner (RHG) model [1]. The bulk modulus of the solid matrix is 5.076 MPsi (35,000 MPa) and that of the fluid is 0.29 MPsi (2000 MPa). Since the problem is a one-unit layer comprising only sandstone, Eq. (28) is applicable. The effective bulk modulus gives a value of 3.7329 MPsi (25,739.05 MPa), which tends towards the bulk modulus of the solid. Figure 2 shows the variation of the bulk modulus with porosity, which indicates an inverse relationship consistent with earlier studies [4,14,19]. Using the random number generator of MS Excel, 50 data points for porosity were created from zero to one and applied to Eq. (28). The trend satisfies the conventional boundary conditions when porosity ranges from zero to one. At porosity equal to one, the value of the effective bulk modulus equals that of the fluid, while at porosity equal to zero, the value equals that of the solid. When the rock is viewed as a unit, the porosity then behaves as an inclusion that reduces the bulk modulus of a solid framework as its concentration increases. Fjær et al. [4] confirmed that the presence of microcracks reduces core samples’ stiffness. These cracks are inclusions into the rocks and so, they reduce the elastic properties of the rock. If the conventional models for deformation take the microcracks into account, it will be extremely hard to define their positions relative to one another. Thus, one cannot overemphasize the advantage of the approach of this study. The porosity equal to zero is analogous to inclusion concentration equal to zero, at which point the rock is strongest. Thus, this study can be useful when engineers aim to quantify the effect of natural fractures (inclusions) on the bulk modulus during drilling waste injection, for example. For such processes, the zones with a large number of natural fractures will possess lower elastic stiffness.

Fig. 1

The wireline log for the Nembe Well-X: Cal is caliper (1 in. = 0.0254 m), Den is bulk density (1 g/cc = 1000 kg/m3), TT is transit time (1 μs/ft = 3.28 μs/m), and GR is gamma ray (API). For depth, 1 ft = 0.3047 m.

View large Download slide

The wireline log for the Nembe Well-X: Cal is caliper (1 in. = 0.0254 m), Den is bulk density (1 g/cc = 1000 kg/m³), TT is transit time (1 μs/ft = 3.28 μs/m), and GR is gamma ray (API). For depth, 1 ft = 0.3047 m.

Fig. 2

View large Download slide

Variation of the bulk modulus with the porosity

The following section compares the model developed in this study with the conventional models for predicting the effective bulk modulus. Using the data in Fig. 2, the models of Voigt, Hashin and Shtrikman, Reuss, and Voigt–Reuss–Hill were used to calculate the bulk modulus. The plots of the bulk modulus against porosity are shown in Fig. 3. The values obtained are between Voigt and Ruess’ models, over the range of the porosity considered, except at the boundaries where all the models converge. The study model (K model) yielded results between the upper and lower bounds and was compared favorably with the HS model. If the extra fluids in the microcracks were not incorporated into the analysis, the study model would yield values similar to Voigt’s upper bound. The distribution of the models under investigation does not follow the same path. This is because of differences in the underlying principles that form their basis. The Voigt model indicates a more linear path than the others as indicated in Fig. 3, while the Ruess’ model dropped significantly because of its tendency towards the bulk modulus of the fluid. If the difference between the modulus of the solid and fluid is narrow, the Ruess’ model will converge towards the other results.

Fig. 3

View large Download slide

Bulk modulus versus porosity for various models

To consider the accuracy of the model of this study, one can use the numerical method to obtain the parameter in Eq. . The equation is recast into the form of a numerical analysis so that one can obtain the value of the parameter (a) that makes the function equal to zero. Then the following expression emerges:

K - (φ + a (1 - φ)^{2 / 3}) K_{f} + (1 - (φ + a (1 - φ)^{2 / 3})) K_{s} = 0

(36)

A numerical scheme requires the use of representative local or regional data for the bulk modulus. However, one can derive a parametric value by taking the average value of the Voigt, V–R–H, and HS models. Figure 4 shows the root (0.16805) when values of the parameter are selected from −1 to 4. This value can be taken as the initial guess of the Newton–Raphson method to get a more representative root. With the data given in this study, the average value of the effective bulk modulus for the chosen models is 3.4257 MPsi (23,621.03 MPa). With the initial guess of 0.16805, the Newton–Raphson method [23] yields a value of 0.168047 after two iterations. The LP Simplex engine of the Excel Solver also gave the same value. Therefore, the expression for the bulk modulus, using a numerical scheme, can also be represented as follows:

\begin{aligned} K = K_{f} (φ + 0.16805 {(1 - φ)}^{2 / 3}) + K_{s} (1 - φ - 0.16805 {(1 - φ)}^{2 / 3}) \end{aligned}

(37)

Fig. 4

View large Download slide

The root of numerical analysis for the Voigt, V-R-H, and HS models

NB: it is more appropriate to use regional or local data to obtain the parameter for any numerical scheme. Figure 5 shows how the study model compares with the numerical model, which was obtained using the average bulk modulus of the selected models.

Fig. 5

View large Download slide

Comparing the model and the numerical bulk modulus

Further validation of the proposed model using FEA follows, but this is not an easy procedure since the study model is not one of the conventional built-in models for such purpose. Moreover, the porosity is not a length-scale for petrophysical analysis [28]; therefore, it cannot be used as a representative elemental volume (REV) for the FEA. Rather, the authors believe that it can be a characteristic of the grids used for the analysis and the solid frame is capable of showing the trend of the bulk modulus. However, the model is valid if it satisfies the linear elasticity of FEA. With a simple 3D geometry characterized by the following: Poisson ratio of 0.3, bulk density of 2348.45 kg/m³, Young’s modulus of the element of 30.88 GPa, and a displacement force of 1000 N, one can validate the model. The force results in work done per volume to displace the poroelastic grids of the frame as follows [5]:

d W = σ_{i j} d ε_{i j} + P_{p} d ξ

(38)

where dW is work done per volume, σ_ij is the total stress tensor, ɛ_ij is the small strain tensor, ξ is a parameter representing fluid content, and P_p is pore pressure. For no change in the fluid content in the grids, the work done reduces to the following expression:

d W = σ_{i j} d ε_{i j}

(39)

Using the externally applied pressure on the element, the work increment per volume for small strain approximates the following expression:

d W = K d ε_{i}

(40)

where the variable ɛ_i is the small volumetric strain in the grid. The integration of Eq. and subsequent application of the pressure-volume work yield the following bulk modulus for the finite elements:

K = \frac{P}{ε_{v}} (1 - φ)

(41)

where P is the compressive stress on the grids, which is obtained from the linear elasticity of the FEATool. Equation is a simple case for the solid frame. Using an REV of 500 m × 500 m × 1850 m for the x, y, and z directions, the authors obtained the following from the computer: 1248 grid cells, which were refined to 2187 grid points with a triangulation of 9987. The grid cell means volume was 549.87 and the grid cell mean the quality of 0.6814. The experiment showed that the displacement of the rock increased with the applied force (Fig. 6), which is consistent with the theory of elasticity [4]. The tool was run under the developed geometry and the applied force was varied to obtain the corresponding displacement. Figure 7 shows the Von Mises stress on the element while Fig. 8 shows the strain caused by the application of force. From the strain in the grids (15.44 × 10¹¹) and the corresponding stress (5.1576 Pa·s), the bulk modulus of the element gave the results shown in Fig. 9.

Fig. 6

View large Download slide

The linear deformation of the elements under bulk compression

Fig. 7

View large Download slide

Von Mises stress on the element

Fig. 8

View large Download slide

The lateral strain on the element

Fig. 9

View large Download slide

Model validation with the FEA

It was stated that the approach of the critical porosity requires slight modification to make it possess the significant characteristics of poroelasticity. In its current state, the model predicts the results of the solid frame although it is capable of yielding results for the porous material as a whole. Equation is different from the conventional relationship between the rock frame and the solid material, which is given as follows [4]:

K_{f r} = K_{s} (1 - φ)

(42)

Although both equations predict that the rock’s frame disappears at porosity equals unity, Eq. defines the absolute characteristics of the solid only, while Eq. implicitly defines the character of the solid and pores. While the model can be modified to include the property of the fluid, the following insight can be added. The grain sizes, sorting, packing, and shape control the porosity [29–31], and as long as the factors causing packing are active, the rock will always be in contact with each other. However, if the sizes of the grains can be made bigger, a point will reach when the grains will escape the dimension of the original frame. This is the maximum porosity, above which the grains will not be in contact with one another inside the original frame. Thus, one can pose the following modified equation that still maintains the requirement of the frame balance:

K_{f r} = K_{s} (1 - {(\frac{φ}{φ_{c}})}^{x})

(43)

Where the parameter x is a correction factor that can be obtained experimentally or numerically. To overcome the challenge of finding the critical porosity, one can apply the boundary condition of the fluid. When the porosity is equivalent to one, the frame approximates the bulk modulus of the fluid; thus, the following expression emerges:

K_{f} = K_{s} (1 - {(\frac{1}{φ_{c}})}^{x})

(44)

Then one can relate the critical porosity to the bulk modulus of the solid and fluid as follows:

{(\frac{1}{φ_{c}})}^{x} = 1 - \frac{K_{f}}{K_{s}}

(45)

The frame bulk modulus now approximates the effective bulk modulus as follows:

K = K_{s} - (K_{s} - K_{f}) φ^{x}

(46)

Note that when the value of the parameter (x) equals unity, one obtains Voigt’s bulk modulus. However, with the use of the data of this study and the numerical scheme, the value of the parameter (x) is 0.8228. Like before, one can apply regional or local data to get a more accurate value for the parameter. Figure 10 is a comparison of the modified frame porosity and the other models. The modified frame model is very close to the average of the HS and VRH models.

Fig. 10

View large Download slide

The display of the modified frame model and the other models

Equations (28), (36), and (46) will most likely yield similar results when the same set of local data is used to obtain the parameters. The Voigt and Reuss models that follow the concept of deformation of layers yield different results most probably because they do not consider the deformation of microcracks. The rate of change of the bulk modulus with porosity depends on the difference between the fluid and solid bulk moduli. This is observed when one differentiates the bulk modulus with respect to porosity. For example, differentiating Eqs. (28), (36), and (46) with respect to porosity shows that the analytical method has the largest recovery of the bulk modulus compared to the numerical and modified frame models.

To validate Eq. , which is the model developed for the crack porosity, use the data from the wellbore and the expression for the secondary porosity. The substitution of Eqs. and into Eq. , and subsequently using the average compaction factor of 1.15 for sonic porosity yield the following:

\begin{aligned} φ_{S} & = \frac{2.65 - 2.348}{2.65 - 1.0} - \frac{73.4 - 55.5}{189 - 55.5} \times \frac{1}{1.15} = \frac{0.302}{1.65} - \frac{17.9}{133.5} \times \frac{1}{1.15} \\ = 0.0664 \end{aligned}

Using Eq. and the computed porosity of 0.163 from the R–H–G equation yield the following:

φ_{c} = (1 - 0.163)^{5 / 3} \times {0.163}^{4 / 3} = 0.7433 \times 0.08904 = 0.0662

The accuracy between these values (over 99%) is an indication that the crack porosity is accurate, yields representative result, and constitutes the major part of the secondary porosity in this wellbore. NB: fracture porosity is commonly lower than 1% of bulk volume in carbonates or clastic rocks [24]. Higher porosity will yield higher crack porosity using Eq. (32).

Conclusion

The goal of this study was to develop a means to predict the bulk modulus of rocks based on the effective medium theory. One objective was to overcome the limitation of conventional models, which do not include the deformation of microcracks. This was achieved by accounting for the volume of fluids in the microcracks and using the effective medium theory. The FEA and petroleum industry equations were used to validate the models developed in this study, and the results of the bulk modulus lie between the bounds of the Voigt and Reuss models. The following points are drawn from the study:

The addition of the fluids in microcracks yields results lower than the Voigt model but higher than the Reuss model.
The presence of microcracks/inclusions reduces the effective bulk modulus of the rock.
When the secondary porosity is mainly from natural cracks, higher primary porosity corresponds to higher crack porosity.
The rate of change of bulk modulus with porosity is higher when one considers the fluids in the microcracks.
The magnitude and concentration of inclusions significantly affect how rocks respond to external stresses.
The critical porosity approach is a sophisticated way of predicting the effective bulk modulus when the boundary condition for fluids is incorporated.

Acknowledgment

Thanks go to Dr. Abrakasa of the Department of Geology, University of Port Harcourt for providing the data from where the authors plotted Fig. 1. Stipends from FUW helped in carrying out the research.

Conflict of Interest

There are no conflicts of interest. This article does not include research in which human participants were involved. Informed consent is not applicable. This article does not include any research in which animal participants were involved.

Data Availability Statement

The authors attest that all data for this study are included in the paper.

References

1.

Nwonodi

,

R. I.

, and

Dosunmu

,

A.

,

2021

, “

Analysis of a Porosity-Based Pore Pressure Model Derived From the Effective Vertical Stress

,”

J. Petrol. Sci. Eng.

,

204

, p.

108727

.

Google Scholar

Crossref

2.

Zimmerman

,

R. W.

,

2017

, “

Pore Volume and Porosity Changes Under Uniaxial Strain Conditions

,”

Trans. Porous Med.

,

119

(

2

), pp.

481

–

498

.

Google Scholar

Crossref

3.

Cheng

,

A. H.-D.

,

2016

,

Poroelasticity

,

Springer

,

Berlin

, p.

16

.

Google Scholar

Crossref

4.

Fjær

,

E.

,

Holt

,

R. M.

,

Horsrud

,

P.

,

Raaen

,

A. M.

, and

Risnes

,

R.

,

2008

,

Petroleum Related Rock Mechanics

, 2nd ed.,

Elsevier

,

Amsterdam

, pp.

219

–

225

.

Google Scholar

Crossref

5.

Detournay

,

E.

, and

Cheng

,

A. H. D.

,

1993

, “Fundamentals of Poroelasticity,”

Comprehensive Rock Engineering

,

J. A.

Hudson

, ed.,

Pergamon Press

,

Oxford

, pp.

113

–

171

.

6.

Berryman

,

J. G.

,

1992

, “

Single-Scattering Approximations for Coefficients in Biot’s Equations of Poroelasticity

,”

J. Acoust. Soc. Am.

,

91

(

2

), pp.

551

–

571

.

Google Scholar

Crossref

7.

Snyder

,

K. A.

,

Garboczi

,

E. J.

, and

Day

,

A. R.

,

1992

, “

The Elastic Moduli of Simple Two-Dimensional Isotropic Composites: Computer Simulation and Effective Medium Theory

,”

J. Appl. Phys.

,

72

(

12

), pp.

5948

–

5955

.

Google Scholar

Crossref

8.

Biot

,

M. A.

,

1941

, “

General Theory of Three-Dimensional Consolidation

,”

J. Appl. Phys.

,

12

(

2

), pp.

155

–

164

.

Google Scholar

Crossref

9.

Wang

,

L.

,

Dresen

,

G.

, and

Rybacki

,

E.

,

2020

, “

Pressure-Dependent Bulk Compressibility of a Porous Granular Material Modeled by Improved Contact Mechanics and Micromechanical Approaches: Effects of Surface Roughness of Grains

,”

Acta Mater.

,

188

(

4

), pp.

259

–

272

.

Google Scholar

Crossref

10.

Fjær

,

E.

,

2009

, “

Static and Dynamic Moduli of a Weak Sandstone

,”

Geophysics

,

74

(

2

), pp.

WA103

–

WA112

.

Google Scholar

Crossref

11.

Fjær

,

E.

,

2019

, “

Relations Between Static and Dynamic Moduli of Sedimentary Rocks

,”

Geophys. Prospect.

,

67

(

1

), pp.

128

–

139

.

Google Scholar

Crossref

12.

David

,

E. C.

,

Fortin

,

J.

,

Schubnel

,

A.

,

Gueguen

,

Y.

, and

Zimmerman

,

R. W.

,

2013

, “

Laboratory Measurements of Low- and High-Frequency Elastic Moduli in Fontainebleau Sandstone

,”

Geophysics

,

78

(

5

), pp.

D369

–

D379

.

Google Scholar

Crossref

13.

Belyaev

,

B. A.

, and

Tyurnev

,

V. V.

,

2018

, “

Electrodynamic Calculation of Effective Electromagnetic Parameters of a Dielectric Medium With Metallic Nanoparticles of a Given Size

,”

J. Exp. Theor. Phys.

,

127

(

4

), pp.

608

–

619

.

Google Scholar

Crossref

14.

David

,

E. C.

, and

Zimmerman

,

R. W.

,

2012

, “

Pore Structure Model for Elastic Wave Velocities in Fluid-Saturated Sandstones

,”

J. Geophys. Res.

,

117

(

B7

), p.

B07210

.

15.

Wang

,

M.

, and

Pan

,

N.

,

2018

, “

Predictions of Effective Physical Properties of Complex Multiphase Materials

,”

Mater. Sci. Eng. R: Rep.

,

63

(

1

), pp.

1

–

30

.

Google Scholar

Crossref

16.

Watt

,

P. G.

,

Davis

,

G. F.

, and

O’Connel

,

R. J.

,

1976

, “

The Elastic Properties of Composite Materials

,”

Rev. Geophys. Space Phys.

,

154

(

4

), pp.

541

–

563

.

Google Scholar

Crossref

17.

Hill

,

R.

,

1963

, “

Elastic Properties of Reinforced Solids: Some Theoretical Principles

,”

J. Mech. Phys. Solids

,

11

(

5

), pp.

357

–

372

.

Google Scholar

Crossref

18.

Zhang

,

L.

,

Ba

,

J.

, and

Fu

,

L.

,

2019

, “

Estimation of Pore Microstructure by Using the Static and Dynamic Moduli

,”

IJRMMS

,

113

, pp.

24

–

30

.

19.

Shapiro

,

S. A.

,

2003

, “

Elastic Piezosensitivity of Porous and Fractured Rocks

,”

Geophysics

,

68

(

2

), pp.

482

–

486

.

Google Scholar

Crossref

20.

Zimmerman

,

R. W.

,

1991

,

Compressibility of Sandstones. Development in Petroleum Science

,

vol. 29

,

Elsevier Science Publishers B.V

,

Amsterdam, The Netherlands

, pp.

1

–

10

.

21.

Hill

,

R.

,

1952

, “

The Elastic Behavior of Crystalline Aggregate

,”

Proc. Phys. Soc. London A

,

65

(

5

), pp.

349

–

354

.

Google Scholar

Crossref

22.

Hashin

,

Z.

, and

Shtrikman

,

S.

,

1963

, “

A Variational Approach to the Elastic Behaviour of Multiphase Materials

,”

J. Mech. Solids

,

11

(

2

), pp.

127

–

140

.

Google Scholar

Crossref

23.

Chapra

,

S. E.

, and

Canale

,

R. P.

,

2007

,

Numerical Methods for Engineers

, 5th ed.,

Tata McGraw-Hill

,

India

, pp.

81

–

118

.

24.

Morton-Thompson

,

D.

, and

Woods

,

A. M.

,

1993

,

Development Geology Reference Manual, AAPG Methods in Exploration Series

. No 10. AAPG, pp.

180

–

209

.

25.

Miskimins

,

J. L.

,

Holditch

,

S. A.

, and

Veatch

,

R. W.

,

2019

,

Hydraulic Fracturing: Fundamentals and Advancements

, SPE Monograph Series, SPE, USA, pp.

1

–

90

.

26.

Oscar

,

G. G.

,

2022

, How to Calculate Fracture Porosity from Well Logs, GeolOil LLC, www.geoloil.com>fracturePorosity, Accessed August 15, 2022.

27.

Vinegar

,

H. J.

,

1986

, “

X-ray, CT, and NMR Imaging of Rocks

,”

J. Pet. Technol.

,

38

(

3

), pp.

257

–

259

.

Google Scholar

Crossref

28.

Knut-Andrea

,

L.

,

2014

,

An Introduction to Reservoir Simulation Using Matlab, User Guild for the Matlab Reservoir Simulation Toolbox (MRST)

,

SINTEF ICT

,

Oslo

, pp.

3

–

10

.

29.

Nwonodi

,

R. I.

,

Ojanomare

,

C.

,

Dosunmu

,

A.

, and

Wami

,

E. N.

,

2020

, “

A Compact and Simple Estimate of Reservoir Rock Tortuosity

,”

Geomechan. Geoeng.

,

17

(

1

), pp.

131

–

140

.

Google Scholar

Crossref

30.

Armitage

,

P. J.

,

Worden

,

R. H.

,

Faulkner

,

D. R.

,

Aplin

,

A. C.

,

Butcher

,

A. R.

, and

IIiffe

,

J.

,

2010

, “

Diagenetic and Sedimentary Controls on Porosity in Lower Carboniferous Fine-Grained Lithologies, Krechba Field, Algeria: A Petrological Study of a Caprock to a Carbon Capture Site

,”

Mar. Pet. Geol.

,

27

(

7

), pp.

1395

–

1410

.

Google Scholar

Crossref

31.

Njoku

,

C.

, and

Pemez

,

C.

,

2011

, Sedimentary Controls on Porosity and Permeability in Deepwater Turbidites, SPE150805-MS, pp.

1

–

5

,

An Equation for the Bulk Modulus of Composites Derived From the Effective Medium Theory

Abstract

Introduction

The Effective Medium for Layered Rocks

Some Conventional Methods Developed for the Bulk Modulus

Methodology

Model Development for the Bulk Modulus of the Composite Rocks

Analysis and Discussions

Conclusion

Acknowledgment

Conflict of Interest

Data Availability Statement

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

An Equation for the Bulk Modulus of Composites Derived From the Effective Medium Theory

Abstract

Introduction

The Effective Medium for Layered Rocks

Some Conventional Methods Developed for the Bulk Modulus

Methodology

Model Development for the Bulk Modulus of the Composite Rocks

Analysis and Discussions

Conclusion

Acknowledgment

Conflict of Interest

Data Availability Statement

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only