The digital rock analysis of biogenically induced reservoir heterogeneities in Cretaceous reservoirs of Saudi Arabia

The characterization of the fluid-flow properties in biogenically altered formations is a key for successful exploration campaigns. This study assessed reservoir quality and evaluated permeability of the biogenically-modified Cretaceous carbonate section in Saudi Arabia. When dealing with these bioturbated carbonates, characterization of sedimentary heterogeneities is often overprinted by the complex spatial geometries of burrows. The high-resolution 3D X-ray microscopy imaging and analysis of these bioturbated sections lead to a better understanding of the interconnectivity between permeable burrows and tight matrix. The analysis deploys multiscale imaging from whole core to sub-micron scale. The coarse imaging at 20-50 microns resolution helped to identify bioturbated sections to select samples for higher resolution tomography. The 3D sample tomograms were segmented to define burrows and matrix distributions, where samples were extracted for thin-sections, scanning electronic microscopy (SEM) and mercury injection capillary pressure (MICP) analyses to refine the pore sizes and rock types. The analysis showed an intricate, highly connected, mixed horizontal and inclined burrow system dominated by Thalassinoides. Intergranular porosity, associated with the fill of Thalassinoides, constitutes a mechanism for permeability enhancement in a tight matrix. Increased permeability is associated with higher dolomite content that might be used as a sweet spot identifier from wireline logs.


Introduction
Textural and mineralogical modifications of the original sedimentary fabric by organisms at the time of sediment deposition are known as bioturbation. Biogenically driven chemical and physical alteration of the primary sedimentary fabric due to burrowing infauna in carbonate formations can result in fabricselective dolomitization [1]. Recent studies have demonstrated that in tight carbonate formations this type of modification can lead to the enhancement of reservoir permeability [1]. Enhancement associated with bioturbation can be expressed as dualpermeability flow-media and shares some conceptual similarities to fracture systems. Sediment modifications in the primary sedimentary fabric include burrowing, churning and continuous reworking of particles. One of the fundamental input parameters into reservoir models is permeability. Commonly, biogenically induced dualpermeability flow media exhibit poor reservoir characteristics owing that only trace fossils (i.e., the permeable conduits) are able to transmit fluids, with little to none interaction with the tighter matrix [1]. This requires detailed understanding of the geological heterogeneities to characterize their impact on petrophysical and reservoir engineering studies. This paper documents the influences of bioturbation on reservoir quality from the Cretaceous Limestone (CL), Saudi Arabia. Permeability distributions were investigated using digital and conventional techniques to refine different rock types. We present a case study of biogenic enhancement in tight facies of a carbonate formation and provide a conceptual framework for burrow-associated permeability improvement analyzed through threedimensional (3D) X-ray computed tomography (CT) imaging acquired at multiple scales for digital rock characterization supported by mercury injection capillary pressure (MICP) experiments and highresolution scanning electron microscopy (SEM) imaging.
The objective of the work was to integrate the CT images, thin-section photomicrographs, MICP and SEM data to simulate permeability in high-resolution 3D images and upscale it to the plug and whole-core levels to define if permeability is linked to specific macroscale geological features of the rock. This may only be examined through multi-scale imaging and modeling. Whole core sample analyses was essential in calculation of the reservoirs heterogeneity to assess permeability anisotropy. The comparison of permeability at different scales is necessary for upscaling laboratory-measured properties to grid cells of the full-field geological models.

Geological context
The CL is a tight formation (<1mD) in Saudi Arabia. To understand the causes of flow restriction, this study captures systematically different pore types, their relationship, distribution, connectivity and their impact on reservoir fluid flow behavior. It is observed that pores are not related to any depositional surface and are rather formed due to mesogenetic corrosion of highly micritized, tight carbonate rock bodies. Primary pores are almost completely destroyed during the process of shallow burial diagenesis. Separated larger pores are both fabric as well as non-fabric selective type. They are the main contributors to pore volume within an overall pervasive (micritic) matrix pore dominated system of wackestone and packstone.

Lithology
Sedimentological analysis based on core samples reveals that burrowed, tight facies range from mudstones to mud-dominated packstones and wackestones with abundant skeletal fragments and pellets.

Depositional environments
Overall, sediment deposition took place on a lowgradient carbonate platform along the eastern border of the Arabian Shield. The preponderance of fine-grained lithology along with intense bioturbation including well-developed Thalassinoides suggests deposition in sheltered lagoonal settings in the inner part of the platform during the Lower Cretaceous.
Several porosity generation and destruction events are recorded in these rocks. A eugenic phase of porosity enhancement took place initially with the dissolution of aragonite shells followed by compaction related porosity reduction. During prolonged subsidence of this passive margin sequence, a good part of the pores were occluded due to calcite cementation.
Pressure solutions are also recorded within argillaceous mudstone, which dissolved carbonate minerals locally and deposited them in the surrounding pores, causing further reduction in porosity. Mesogenic porosity enhancement is found as an important event for this formation. It is presumed that a phase of organic acid generation and carbonate mineral dissolution took place during the onset of hydrocarbon expulsion from deeper source beds. This process was very important for the CL as most of the porosities are developed through secondary pores.

Conventional core observations
Core descriptions, thin-section petrography and SEM assessments allowed the identification of an intricate, highly connected, mixed horizontal and inclined burrow system dominated by Thalassinoides. Selective dolomitization occurs in burrow fills and is commonly surrounded by a non-dolomitized lime mud matrix. This selective burrow dolomitization results in intergranular porosity and constitutes a mechanism for vertical permeability enhancement in an otherwise lowpermeability matrix. Preferential vertical burrowing direction suggests that dwellings of bioturbating infauna constitute selective fluid-flow networks. Core pictures displaying oil-stained Thalassinoides burrows show that oil staining is concentrated preferentially in the fill of burrows ( fig. 1, a). Thin-section photomicrographs provide information about rock texture and porosity types in carbonates that are essential for rock-type classification. Photomicrographs ( fig. 1, b, c) illustrate dolomitized burrows, muddy matrix and intergranular porosity concentrated in the burrow fill due to preferential dolomitization associated to Thalassinoides. Burrows are mainly vertical with a minor horizontal component, and may span a depth of up to one meter below the sea floor, with typical diameters between five and ten millimeters.
Bioclastic material is partially leached and cemented by drusy calcite. This section of the sample shows good porosity with common fine moldic pores and micropores. A calcite cemented fracture separates the dolomitized and non-dolomitized parts of the sample. Microporosity is commonly observed in the bioclastic wackestone.

Scanning electron microscopy
Samples containing burrows and matrix were selected from each whole core sample. Continuous magnification of SEM observations together with energy dispersive spectroscopy (EDS) show that most of the pores preserved in CL is moldic. Images ( fig. 2, a-d) highlight the dominance of solution pores within the formation intervals. Dolomitization of grains selectively increases intercrystalline porosity locally, but this phenomenon is less common. Reservoir quality as a whole is poor due to poor connectivity of secondary dissolution or intercrystalline pores. Burrows have large dolomite grains. Matrix has tight micritic calcite as well as a few large dolomite grains.
The main goal of the SEM-EDS measurements was to establish the presence of dolomite grains in the micritic system. Fig. 2 (e, f) shows SEM-EDS points. The results demonstrate that grains are mainly calcite and dolomite.

Digital rock analysis of the CL bioturbated samples
The digital rock analysis workflow includes the whole core samples computed tomography analysis, burrow segmentation, identification of burrow-modified regions, and the selection of smaller samples for higher resolution work. The smaller samples were scanned at medium resolution (4 μm) to identify the matrix to burrow interface, and at high-resolution (0.5 μm) to study morphology and properties of burrows and matrix. Based on initial image analysis results, appropriate volumes were taken for MICP. Once MICP results were available, sub-samples that cover the entire porosity spectrum in the rock were selected and extracted. The next phase of the workflow involved the determination of burrow and matrix properties through high-resolution multi-scale CT imaging. In addition, a combination of high-resolution CT and SEM-EDS data were used to refine the mineralogy of the bioturbated zones. Image processing and various analyses lead to computations of porosity and permeability of the burrowed and matrix sections of the rocks. The final stage of the workflow is to integrate the properties of burrows and matrix to produce up-scaled properties. These results can also be integrated to the whole core and wireline log scale; however, this is beyond the scope of this paper.

Whole core samples imaging and analysis
The whole core samples were selected based on the sedimentological descriptions of 200 ft of the whole core and 1,000 ft of well logging data. Sampling locations were identified to represent variations in lithofacies, texture and porosity type. In total, thirty 25 cm long whole core samples were taken for lowresolution (25 μm) scanning. Examples of the whole core samples scans are illustrated in fig. 3.
The grayscale image of the whole core samples ( fig. 3) shows that it is composed of two distinct lithology regions: matrix and longitudal circular shaped Thalassinoides burrows.
Based on the low-resolution scanning three type of sub-samples have been selected for the following studies: advanced high-resolution computed tomography (16 samples), MICP measurements (8 samples) and conventional core analysis to define total porosity and absolute permeability (37 samples). Visual examination of all samples was carried out to ensure the plugs were physically suitable for the measurements and free of visible fractures or anomalies.
Sub-samples were cut from representative regions in the core that preserved changes in the geological features seen in the whole core CT images to provide integrated characterization among available data. This is an effective way of understanding rock structure at multiple scales, which would lead to proper evaluation of the impact of rock quality and geological heterogeneity on the permeability of the reservoir rocks.

Mercury injection capillary pressure data analysis
The laboratory-measured MICP samples are carefully cut from the whole core to represent each phases (burrow and matrix) in the sample (usually 1 inch diameter plug).
The MICP data indicate the prevalence of micropores in matrix ( fig. 4, a) and mesopores in bioturbated regions ( fig. 4, b) of the formation. MICP data show a narrow, unimodal distribution of the porethroat radii, highlighting the development of mesopores in burrows, with the moderate entry pressure and gentle saturation gradient indicating fluid flow via mesopores. MICP data observations on the matrix sections indicate that pore-throat size distribution is bimodal, representing two types of micropores: 0.03 μm and 0.3 μm in diameter.
The difference in the pore-throat diameters is interpreted as the cause of higher permeability in the burrowed sections. Absolute permeability variation associated with difference in the pore-throat size diameters has been investigated with conventional core analysis (CCAL) of the plugs.

Conventional core analysis of plugs
Plug samples (1.5 inch diameter) were taken from each lithology to represent the statistical property variations along the entire core length. The samples were cleaned in hot solvent extractors and were measured for helium porosity and gas permeability in standard equipment at ambient conditions using a calibrated steady-state permeameter with nitrogen gas as flowing medium.
CCAL demonstrates clear separation between bioturbated and matrix intervals on porositypermeability graph (fig. 5). Matrix is tighter and exhibits a smaller range of porosities in comparison to the burrows porosity range. The porosity in this burrowed phase is the main contributor to the flow. On average, total porosity of these tight bioturbated plugs exceeds matrix porosity by 5-7% and permeability difference consists 1-2.5 orders of magnitude.
CCAL results underscore the importance of the bioturbated sections analyses; however, it does not provide information on permeability (k) anisotropy in vertical (k v ) and horizontal (k h ) directions. Another important task was to assess the impact of the burrow to matrix volume relationship with bulk permeability of the sample. High-resolution digital rock analysis was utilized to answer these questions.

High-resolution digital rock analysis
Multi-scale digital rock analysis is conceptually wellsuited for bioturbated rocks analysis and modeling [2] as it allows varying resolution and provides a systematic procedure for coarsening and refinement. Darcy's model is to approximate pressure and fluxes on a coarse grid in large-scale discontinuities, whereas fine-scale effects are captured through basis functions computed numerically by solving local Stokes-Brinkman flow problems on the underlying fine-scale geocellular grid. The Stokes-Brinkman equations give a unified approach to simulating free-flow and porous regions using a single system of equations, they avoid explicit interface modeling, and reduce to Darcy or Stokes flow in certain parameter limits [3]. The modeling demonstrates how fine-scale flow in bioturbated pore networks can be represented within a coarse-scale Darcy-flow model in fractures of the burrow to matrix interface and multi-scale elements. The flow computed solving the Stokes-Brinkman equations in burrow and matrix. This method is a promising path toward direct simulation of highly detailed rock models in bioturbated formations. The high-resolution CT imaging was conducted to build pore structures for the flow analysis.

Imaging and segmentation of full-plugs
To perform high-resolution CT imaging, 1.5 inch samples that include a combination of burrow and matrix were extracted and scanned at 4 μm resolution. Orthogonal slices of burrow and matrix at this resolution are shown in fig. 6. It illustrates that small pores are not present at this resolution ( fig. 7) within the matrix. The throat sizes distribution derived from MICP data (see fig. 4) shows that for determining the flow properties of these burrow-modified cores, the required resolution is at least 5 μm for burrow-throughflow and below 1 μm for the matrix flow unit. The higher resolutions scans are typically acquired to understand the fabric of the rock, and to define largescale geologic objects in the samples. Therefore, a range of resolutions is needed to fully describe the plugs.

Imaging and segmentation of sub-plugs
The analysis of the burrow flow unit requires drilling of sub-plugs in the matrix to burrow interface locations and scanning them at 1 μm resolution; smaller separate burrow volumes were analyzed at maximum 0.5 μm that revealed a connected pore system in the burrow phase. Figure 8 shows a cross section CT slice for the burrow to matrix flow unit. The digital rock analysis confirmed the SEM data. The large dolomitic grains can be seen predominantly in the burrow phase, and the space between them is filled with micritic calcite ( fig.  9). Based on the MICP data, the matrix flow unit would require imagining by focused ion beam SEM beyond the resolution of CT images. Therefore, permeability of matrix for pore-scale modeling was assigned based on the minimum CCAL of samples with matrix only. The modeling results will represent the pessimistic estimate of the whole system permeability.

Numerical modeling of the flow properties
Numerical modeling of permeability in burrows was based on the highest magnification tomography at 0.5 μm resolution. 3D images were segmented to allow porosity and permeability computation at subsample scales. Two main phases have been identified: pores and grains.
Incompressible flow in a porous rock matrix typically obeys Darcy's law and is described by a firstorder elliptic system in which Darcy's law was combined with a mass-conservation equation to relate the pressure and the total (interstitial) velocity.
Incompressible flow in open domains, on the other hand, obeys the Stokes equations [3,4]. The Stokes-Brinkman equations combine Darcy and Stokes models into a single equation. This model gives a unified approach to model flow in both the fractures and the intergranular porous subdomains using a single system of equations. In the free-flow (or fluid) domain, it is assumed that permeability tends to infinity and sets the effective viscosity equal to the fluid viscosity, otherwise it transforms into the coupled Darcy-Stokes equations, which reintroduces the requirement for interface conditions and computational intractability.
The original structure, porosity segmentation, velocity field and illuminated streamlines are illustrated in fig. 10. The porosities of burrows are in the range of 14 to 19 pu, while the permeability is in the range of 5.25 to 2.55 mD for the studied samples that compares well with literature sources [5].

Modeling of the fluid flow through composite rocks
Flow-based upscaling procedure based on Stokes-Brinkman equation computes effective permeability on coarser grids. To build the coarser grid, it is required to estimate permeability of the unresolved porous materials within the grid. Three sub-samples were selected to characterize the burrow and matrix system under sub-micron CT imaging ( fig. 11): a) burrow phase, b) burrow to matrix interface and c) matrix.
As discussed before, permeability of 0.01 mD was pessimistically assigned to matrix based on CCAL results of the typical matrix samples. Burrow permeability was estimated with Stokes-Brinkman equation in 2.4.3.
The extracted properties of burrow and matrix then were combined into the larger scale network to evaluate permeability and interconnectivity. Two composite models were tested: 1) direct communication between burrows and matrix, and 2) inter-burrow flow through the fracture interface. In both cases, the fluid-flow was modeled with Stokes-Brinkman numerical solver under no-slip boundary conditions.
Unresolved porous burrow and tight matrix fabrics have been propagated in 3D volumes through utilization of the machine learning algorithm based on the Extra Trees engine [6] complemented with resolved porosity. This classifier implements a meta estimator that fits a number of randomized decision trees on subsamples of the burrow and matrix grayscale dataset and uses averaging to improve the predictive accuracy and control over-fitting. Examples of the original grayscale image and classification results are shown on fig. 11 and fig. 12.
The results of the modeling illustrate permeability enhancement and porosity-permeability relations in the bioturbated system in three orthogonal directions. The interaction of the burrow system with the fracture network contributes to the bulk enhancement of flow between dead-end burrows and improves permeability of the system by 1-2 orders of magnitude in comparison with the case of direct communication between burrow and matrix ( fig. 13).   CCAL confirms enhancement of the system permeability by the fractures in the burrow to matrix interface. Highest permeability CCAL samples have been scanned and demonstrate various degrees of the fracturing. The porosity-permeability graph in fig. 14 shows five fractured samples together with the 3D render of one of them.
Comparison of the CCAL data against digital rock analysis results (see fig. 14) shows general match between two datasets with commonly underestimated porosity of matrix samples due to limitations of the scanning resolution not showing the microporosity smaller than 0.5 μm.

Relation between burrow volume fraction and permeability
The developed digital rock model was used to simulate the permeability ratios K v /K h as a function of burrow and dolomite volume fractions shown on fig. 15. The volume fractions of burrows vary from 0.2 to 0.75 and volume fractions of dolomite vary from 0.15 to 0.65 that is comparable with [7]. Increase of both burrow and dolomite fractions content is leading to higher K v /K h ratios. These relations ( fig. 15) can be used for the permeability anisotropy prediction based on dolomite content derived from core or well log data analysis to improve permeability simulation in the grid of geological models.

Conclusions
Multiscale digital rock analysis was completed on whole core samples, full plugs and high-resolution subplugs to investigate the geological and petrophysical properties of the bioturbated CL at varying resolution from 25 μm to 0.5 μm. Integrated analysis of the CT images, thin-sections, CCAL and MICP data resolved permeability of burrow and matrix as well as permeability ratio in vertical and horizontal directions.
In the case study presented herein, high burrow connectivity suggests that a continuous flow network is likely to be established. The analysis shows, this is primarily influenced by the: 1) concentration and preferential orientation of permeable trace fossils; 2) degree of permeability contrast between burrows and surrounding matrix; and 3) three-dimensional connectivity of the Thalassinoides network system. Here, Thalassinoides occupy a significant volume of the rock compared to the fractures and this fact enhances the permeability and transmissivity of otherwise low-permeability matrix.
Detailed thin-section descriptions and MICP experiments improved the understanding of the geological variations in the 3D rock volume. This analysis has led to an improved geological interpretation of the diagenetically altered carbonate rocks. The overall analysis offered improved understanding of the heterogeneity in the formation and its impact on multi-scale permeability measurements.
There is at least one order of magnitude difference in permeability of burrow and matrix. Permeability of matrix appears isotropic. Permeability of the burrow fraction is one to three orders of magnitude larger than the matrix permeability; the results show enhanced burrow permeability with higher dolomite fractions. The flow between dead-end burrows occurs through the matrix or fracture interface as presented by the composite models. The K v /K h ratio as a function of burrow and dolomite volume fractions are identified. The ratio of K v /K h is always above one for bioturbated samples.
Increased permeability associated with the higher dolomite content resulting from intensely bioturbated intervals might be used as sweet spots identifier from wireline logs. Recognition of bioturbated intervals and the role they play in reservoir petrophysics may aid in both the optimization of exploration and field development campaigns in similar biogenically modified intervals in the study area.