## 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

*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. (1) to represent the volumes of each layer as follows:

*V*

_{1}is the volume of layer 1,

*V*

_{2}is the volume of layer 2,

*h*

_{1}is the thickness of layer 1, and

*h*

_{2}is the thickness of layer 2. The total volume of the layers comprises the sum of the contributions from each layer:

*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:

*K*

_{m}is the effective bulk modulus,

*K*

_{1}is the bulk modulus of component 1, and

*K*

_{2}is the bulk modulus of component 2.

## Some Conventional Methods Developed for the Bulk Modulus

*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. (7), 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.

*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}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*

_{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. (11) to finally arrive at Eq. (9), assuming the fluid occupies the pores only. From Eq. (11), the bulk compressibility will tend to a uniform value if the components have approximately equal compressibility.

*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.

*φ*

_{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.

*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

*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*

_{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

*K*

_{f}is the bulk modulus of the fluid component, $%Vs$ 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. (17). 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}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:

*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.

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.

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

*K*

_{fi}is the bulk modulus of the fluid in the

*i*th 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. (29).

*φ*

_{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]:

*ρ*

_{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.

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.

*a*) that makes the function equal to zero. Then the following expression emerges:

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.

^{3}, 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]:

*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:

*ɛ*

_{i}is the small volumetric strain in the grid. The integration of Eq. (40) and subsequent application of the pressure-volume work yield the following bulk modulus for the finite elements:

*P*is the compressive stress on the grids, which is obtained from the linear elasticity of the FEATool. Equation (41) 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

^{11}) and the corresponding stress (5.1576 Pa·s), the bulk modulus of the element gave the results shown in Fig. 9.

*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:

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.

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.

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.