Printer Friendly

Modeling of multiscale porous media.


Reservoir rocks are highly heterogeneous and many of them contain a complex pore structure with a wide range of length scales. One of the main objectives of pore scale physics research is to predict macroscopic petrophysical properties from the underlying porous microstructure. Apart from the basic understanding of the physics of flow processes, this has also enormous practical importance in terms of improving uncertain estimates of recovery efficiency on larger scales. Three dimensional pore scale models of the rock are essential for such research (Adler, 1992; Hilfer, 2000; Thovert et al., 2001; Arns et al., 2003; Okabe and Blunt, 2007).

Among the frequently occurring multiscale porous media in nature, carbonate sediments are of great economical interest because nearly half of the worlds petroleum is found in carbonate reservoirs. Carbonate rocks in general, and dolomites in particular, are formed by a large number of different physicochemical processes and contain a wide variety of diagenetic textures overprinting their primary facies. Frequently, the original primordial depositional textures (e.g., calcite ooids or bioclasts etc.) remain visible in the final microstructure (Moore, 1989; Lucia, 1999; Moore, 2001). They contain vuggy pores, defined as pores that are connected only through interparticle porosity (Lucia, 1999). Dissolution, nucleation and growth of crystallites are common phenomena during dolomitization processes, and both porosity and texture can fluctuate in the evolution of the rocks. Unlike sandstones, permeability in carbonate rocks can vary widely within a given rock sample (Fernandes et al., 1996; Anselmetti et al., 1998). Due to these reasons it has been difficult to fully classify and characterize the pore scale microstructure of carbonate rocks (Dunham, 1962; Lucia, 1999).

The geometrical and petrophysical parameters of carbonates depend strongly on resolution because of a very wide range of pore and crystallite sizes dominating its geometric texture. For typical carbonate rocks, with smallest calcite or dolomite crystallite sizes of [10.sup.-7] m, one needs at least a voxel size of a [approximately equal to] [10.sup.-8] m to resolve them. Even for a small cubic sample of size L = [10.sup.-2] m, the digitized sample requires prohibitive CPU time and storage space (~ [10.sup.18] digits). To be able to model such a system and at the same time capture the complex geometry, a completely different modeling approach and data structure are needed. We propose here a pore scale model in the continuum that overcomes these impediments and can act as a starting point for more realistic modeling of these multiscale porous media.

Because the diagenetic processes that produce the large variety of carbonate rock textures are complex and largely unknown, we attempt a simple stochastic geometrical model that tries to reproduce the main features of the pore scale geometry. Any modeling method for carbonates requires solving two key problems: the "pore sizes" and grain diameters ranging over many decades in length scales and the carbonate textures showing a strongly correlated disorder (Song and Sen, 2000). It is also crucial to correlate pore scale microstructure to the depositional and diagenetic fabrics and simultaneously distinguish between intergrain, intragrain, intercrystalline and vuggy pores. Porosity depends strongly on the shape and size of the grains and grain packing. The sizes of intergrain pores depend on the grain size sorting and decrease as cement occludes the pore space or as grains are forced closer together by compaction (Lucia, 1999).

The model proposed here aims to incorporate simultaneously the above diagenetic processes and the following microstructure features: scale dependent intergranular porosity over many decades, vuggy porosity, a percolating pore space, a fully-connected matrix space, intergranular correlations from primordial depositional textures and mineral transformation. It allows discretizations at any arbitrary resolution. In this paper the feasibility of this model is tested by reconstructing the generic microstructure of two different carbonate rocks with multiple length scales.


In this continuum growth model the porous medium is represented by a random sequence of points decorated with objects (grains or crystallites) that correlate with the location and properties of the different crystallite phases present in the original rock. This model belongs to the class of so-called germ-grain models (Stoyan et al., 1995). The deposition of the points is similar to a random sequential adsorption process (Feder, 1980). A stochastic geometer might call the model a random-field-controlled germ-grain model of random sequential adsorption type.


The state space of a rock sample containing N crystallites occupying a bounded region S[subset] [R.sup.3] is the set

[[OMEGA].sub.N] = [(S x [[R.sub.min], [R.sub.max]] x E x {1, 2, ..., g}).sup.N] (1)

of all sequences [omega] = ([[omega].sub.0], [[omega].sub.1], [[omega].sub.N], [member of] [[OMEGA].sub.N]) with [[R.sub.min],[R.sub.max]] [subset] [R.sup.1] and E = {x [member of] [R.sup.3] : [approximately equal to] = 1}. The sample contains g different crystallite phases. An element [[omega].sub.i] = ([x.sub.i],[R.sub.i], [a.sub.i], [T.sub.i]) of the sequence represents a crystallite of type [T.sub.i] at spatial position [x.sub.i] [member of] S with a radius of the equivalent sphere [R.sub.i] and orientation [a.sub.i]. These crystallite attributes shall depend on the original depositional texture through their distribution. A probability distribution Pr on the space [OMEGA] of sequences further specifies the model.


The crystallite phases and intergranular correlations from the underlying primordial depositional textures are first captured by a greyscale image G defined mathematically as a bounded, but not necessarily continuous function G : S[right arrow] [0,1]. This crucial input function G(x) is constructed with information from geological analysis of the rock, quantitative image analysis of two dimensional micrographs and three dimensional [mu]-CT images, if available.

The location and properties of the crystallites in [??] are then defined through:

[R.sub.i] = R(G([x.sub.i])), [a.sub.i] = A(G([x.sub.i])), [T.sub.i] = T (G([x.sub.i])), (2)

where the functions R : [0, 1] [right arrow] [[R.sub.min],[R.sub.max]] [union] {0}, A [0,1] [right arrow] E and T : [0,1] [right arrow] {1,2,...,g} specify their correlations with the primordial texture. The permitted crystallite sizes in the interval [[R.sub.min],[R.sub.max]] correlate with the polydispersity in each crystallite phase of the original sample. The rock sample is represented as a random sequence of N points, each decorated with crystallites satisfying the primordial correlation through Eq. 2.


The points [x.sub.i] are chosen randomly in S C [R.sup.3] and added sequentially if the following overlap condition is satisfied. For each chosen point [x.sub.i] a corresponding sphere radius [R.sub.i] is chosen as defined in Eq. 2 and we set Pr([[OMEGA].sub.c],) = 0 for

[[OMEGA].sub.c] = {[omega][member of][OMEGA]:[there exists]i,j,O([[omega].sub.i],[[omega].sub.j])[??](0,[[lambda].sub.i]]} (3)

where the degree of overlap between the spheres with radii [R.sub.i] and [R.sub.j] is measured by

O([[omega].sub.i],[[omega].sub.j]) = [R.sub.i] + [R.sub.i] - |[x.sub.i] - [x.sub.j]|/[R.sub.i] + [R.sub.j] - |[R.sub.i] - [R.sub.j], (4)

and O([[omega].sub.i], [[omega].sub.j]) is set to zero when [R.sub.i] or [R.sub.j] or both vanish. The parameter [[lamda].sub.i] = [LAMBDA](G([x.sub.i])) defines the maximum allowed overlap with the crystallite at [x.sub.i] and [LAMBDA] : [0,1] [right arrow] [[[lambda].sub.min],[[lambda].sub.max]] with 0 < [[lambda].sub.min] [[lambda].sub.max] < 1. In other words, each newly added crystallite has a finite overlap with an existing crystallite where the degree of overlap is defined by the primordial filter function. This ensures full matrix connectivity. Porosity and pore space connectivity of each crystallite phase depends on the degree of overlap, the density of crystallites and the corresponding crystallite size distribution. The density of crystallites depends on the density of points, denoted by [rho], deposited in the corresponding phase. This is a model parameter that is determined from the image analysis of the original rock. Porosity and pore space connectivity of the full sample depends on the original pore phase and vugs defined in G, inter-crystallite pore space in each phase and the way these different crystallite phases combine to form the full rock sample.


Apart from inter-grain and inter-crystalline pore space, vuggy pore space is a distinctive feature of many multiscale porous media. In general, vuggy pores may be touching or non-touching and connected by other type of porosity (Lucia, 1999). To include such vuggy pores, deposition of crystallites at vuggy pore regions are restricted by setting Pr([[OMEGA].sub.v]) = 0, for

[[OMEGA].sub.v] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N; R(G([x.sub.i]))=0} (5)

This correlates the vugginess with the depositional texture.


After depositing random points in a correlated fashion, the deposited points are then decorated with crystallites according to the depositional texture (Eq. 2). A crystallite of type [T.sub.i] with orientation [a.sub.i] is placed at each point [x.sub.i]. For a feasible and fast convergence of the deposition algorithm, the overlap condition in Eq. 3 is defined for spherical volumes associated with each deposited point. Since the decorated crystallites can be of arbitrary shape, the spherical volume must either be inscribed within the crystallite or the size of the crystallite [d.sub.i]([R.sub.i]) must be scaled appropriately so as to retain the matrix connectivity.


The porous media sample is fully represented by a list of N quadruples ([x.sub.i],[d.sub.i],[a.sub.i],[T.sub.i]). It can be discretized at any arbitrary resolution. In the discretization procedure, a cubic sample of sidelength l is subdivided into a grid of cubic voxels, each of side length a. Within each voxel a set of selected collocation points are chosen and the voxel is then assigned a value depending on the fraction of these collocation points that fall inside the deposited crystallites. A larger set of collocation points increases the accuracy of the discretization, but is computationally demanding when N is large.


In this simple but reasonably effective discretization procedure a voxel is marked as matrix and assigned the label i if all of the following nine points fall within the ith crystallite [G.sub.i], i = 1, 2,...,N.

[p.sub.j]=p+a/2([e.sub.1]+[e.sub.2]+[e.sub.3])+a/4[t.sub.j]; j=0,1,..,8, (6)

where p is the position of the first corner of the voxel, a is the sidelength of the voxel, [t.sub.0] = (0,0,0), [t.sub.1] = (1,1,1), [t.sub.2]=(-1,1,1), [t.sub.3] = (1,-1,1), [t.sub.4] = (1,1,-1), [t.sub.5] = (-1, -1,1), [t.sub.6] = (1,-1,-1), [t.sub.7] = (-1 ,1,-1) and [t.sub.8] = (-1,-1,-1). These nine points are the center and the eight symmetric points between the center and the corners of the voxel. So, if the voxel is fully inside a single crystallite then its label is the crystallite number. If all the nine points fall inside the pore space, i. e., outside all the N crystallites, then the voxel is resolved as pore and assigned a value 0. In all other cases the voxel status is labeled as unresolved at the current resolution and assigned the label -m if m points fall inside one or more crystallites. This label reflects the matrix density within the voxel. Such labeling criteria with crystallite information is useful for building network models and computing mechanical properties from the discretized representation. The above discretization rule with just nine collocation points provides a reasonably good discretization for the analysis of the model. A more accurate discretization method resembling experimental [mu]-CT procedure is presented below.


The number of collocation points in the discretization rule is increased from 9 to [n.sup.3]. These [n.sup.3] collocation points are the points of a n x n x n cubic sublattice positioned centrally inside each voxel. Each voxel is labeled with an integer number m with the following meaning

- m = [n.sup.3] : all collocation points are in matrix space. - m = 0 : all collocation points are in the pore space. - 0 < m < [n.sup.3] : exactly m collocation points are in the matrix space.

This procedure produces synthetic [mu]-CT discretizations of the computer model of the rock. Experimental [mu]-CT images contain noise and a direct comparison may require addition of noise to the synthetic [mu]-CT or filtering of the experimental [mu]-CT or both. Digitized representations of the pore scale can be obtained by choosing an appropriate threshold [m.sub.c] and relabeling the voxels with label 0 < m < [m.sub.c] to 0 (pore) and the voxels with label m, [less than or equal to] m [less than or equal to] [n.sup.3] to 1 (matrix). It is important to note that the CPU time needed for the discretization procedure is proportional to the number of collocation points ([n.sup.3]) and can be prohibitive if the number of crystallites defining the microstructure is large.


The discretization procedure requires us to check if a given point p [member of] [G.sub.i], i = 1, 2,...,N. The computational procedure for dolomite type crystallites, i. e., [G.sub.i] is rhombohedral, is as follows: A dolomite crystallite [G.sub.i] at position [x.sub.i] is defined by three pairs of intersecting parallel planes. Each pair is separated by a distance [d.sub.i] and tilted by an angle [alpha] degrees about the coordinate axes to which they are parallel initially. The equations of these three pairs of planes are

[n.sub.j].x[+ or -][d.sub.i]/2=0, j=1,2,3, (7)

where [n.sub.1]= (cos [alpha], 0, - sin [alpha]), [n.sub.2] = (- sin [alpha], cos [alpha], 0), [n.sub.3] = (0,-sin [alpha], cos [alpha]).

The rotations of the crystallites are implemented through generalized quaternions. For each [G.sub.i], with the center of the crystallite at the origin of the coordinate system, the orientation [a.sub.i] is defined by a sequence of three rotations of [[theta].sub.1], [[theta].sub.2] and [[theta].sub.3] about the coordinate axes [e.sub.1] = (1,0,0), [e.sub.2]= (0,1,0) and [e.sub.3] = (0,0,1) respectively. So

[a.sub.i] = A(G([x.sub.i])) = -[Q.sub.i3][Q.sub.i2][Q.sub.i1], (8)

where the unit quaternion [Q.sub.ij] represents a rotation of [[theta].sub.j] about the vector [e.sub.j]

The point p [member of] [G.sub.i] if

([n.sub.j].[p'.sub.r] + [d.sub.i]/2)([n.sub.j].[p'.sub.r]-[d.sub.i]/2) < 0, j=1,2,3, (9)

where [p'.sub.r] = [q.sub.i.sup.-1] (p-[x.sub.i])([q.sub.i.sup.-1])* and (.)* is quaternion conjugation. The case [alpha] = 0 corresponds to cubic [G.sub.i]. Checking p [member of] [G.sub.i] for spherical shaped grains is straightforward. The above procedure can be readily generalized to other arbitrary shaped crystallites defined by cutting planes.


Realistic models of multiscale porous media can be generated using the algorithm described in Eqs. 1-5. Actual implementations may differ from sample to sample depending on the underlying primordial textures, crystallite size distributions, type of crystallites, location and type of vugs etc. Apart from testing the feasibility of the microscopic modeling for carbonate rocks, the following description presents the detailed computational implementation of this modeling procedure for two different carbonate microstructures.



Here we present in detail the specific reconstruction of an oolithic dolostone (Biswal et al., 2007). Only a two-dimensional photomicrograph of the original sample was available as input data in this case. A part of the image is shown in Fig. 1. It contains ellipsoidal primordial facies, complete dolomitization and vuggy porosity of touching-vug type.

Image analysis

Due to insufficient resolution and small size of the original image, computerized quantitative image analysis is not feasible. Therefore, only visual inspection and manual measurements are carried out. The rock contains ellipsoidal ooid grains that are fully dolomitized to two distinct crystallite phases: a thin isopachous layer (rims of the elliptical ooids) and the inner core. The diameters of the dolomite crystallites in these two phases are roughly in the interval (5 [micro]m to 25 [micro]m) and (1 [micro]m to 10 [micro]m) respectively. No discernible crystallite size distribution or overlap variation are visible in both phases. The crystallite packing in both phases shows a space filling by large number of smaller crystallites around bigger crystallites. The crystallite packing in the inner cores of ooids have large overlap. The ooid grains feature diagenetic replacive dolomitization and dissolution vugs. Thresholding analysis of the image indicates the porosity of the sample to lie in the range of 0.25-0.3.

Primordial filter function

In the absence of any three dimensional information, the primordial filter function G(x) is constructed from a computer model of polydisperse sedimentation of ellipsoidal ooids (disc-shaped) in S [subset] [R.sup.3] (Bakke and Oren, 1997). Vuggy pores were artificially generated by removing a small fraction of the deposited ooids and enlarging the remaining ooids by 20% of their size to compensate the inter-ellipsoid pore space and porosity. The kth ooid represented by {[r.sub.k],[s.sub.k],[e.sub.k]} is centered at [r.sub.k] = ([x.sub.k],[y.sub.k],[z.sub.k]) with semi-axes lengths [a.sub.k] = [e.sub.k][s.sub.k], [b.sub.k] = [e.sub.k][s.sub.k], [c.sub.k] = [s.sub.k] and k = 1, 2, ..., 1359. The primordial filter function is defined as


where [d.sub.k] = [[[([x.sub.i]-[x.sub.k]).sup.2]/[e.sup.2.sub.k] + [([x.sup.2]-[y.sub.k].sup.2]/[e.sup.2.sub.k] + [([x.sub.3]-[z.sub.k]).sup.2]/1].sup.1/2]. A small thin section of the primordial geometry is shown in Fig. 2 (top left).


Deposition of points

The reconstructed cubic carbonate rock sample for Sample 1 has a sidelength l = 2 mm, i. e., S = [[0,l].sup.3]. Due to large polydispersity and [R.sub.min]/l << 1, the continuum list contains millions of crystallites. The deposition rule requires a computational scheme for random deposition of overlapping polydisperse spheres in S that obey a specific size distribution and each added point [x.sub.i] [member of] S must satisfy Eq. 3.

For computational efficiency we divide S [subset] [R.sup.3] into smaller non-overlapping cubic cells U each of sidelength l, i. e., S = [U.sub.1 [U.sub.2] [U.sub.3]...Points are deposited randomly in each cell Ui and added to the final list if Eq. 3 is satisfied. This scheme provides fast convergence in deposition as Eq. 3 is verified only for a small set of existing points in [U.sub.i] and the neighboring cells. Further, to ensure a space filling by a large number of smaller crystallites around bigger crystallites as observed in the original image, the following crystallite size distribution is incorporated into the deposit scheme. For the jth deposited point at [x.sub.j] in each cell U, the excluded volume radius [R.sub.j] is chosen as


where [DELTA]R = [R.sub.max] - [R.sub.min] and [[xi].sub.j] is a uniform random number in (0,1). Although [[lambda].sub.i] = [LAMBDA](G([x.sub.i])), for this sample, we have assumed it to be constant, [[lambda].sub.i] = [[lambda].sub.c],. In each cell n = [rho][l.sup.3] number of points are deposited. The packing density p is first determined from trial samples that match the target porosity and show pore space connectivity. This is checked using the HoshenKopelman algorithm (Stauffer and Aharony, 1992) on high resolution discretizations of the trial samples.

For the crystallites in the isopachous layers, points are added to S with Pr([[OMEGA].sub.1]) = 0, for

[[OMEGA].sub.1] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, G([x.sub.i]) [??] (0,0.125)} (12)

with [R.sub.min] = 2.5 [micro]m, [R.sub.max] = 12.5 [micro]m, [[lambda].sub.c] = 0.6, l = 62.5 [micro]m and [rho] = 1.95 x [10.sup.-4]. For crystallites in the inner core of the ooid grains points are added to S with Pr([[OMEGA].sub.2]) = 0, for

[[OMEGA].sub.2] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, G([x.sub.i]) [member of] [0,0.1]} (13)

with [R.sub.min] = 0.5 [micro]m, [R.sub.max] = 5 [micro]m, [[lambda].sub.c] = 0.65, l = 7.5 [micro]m and [rho] = 3.395 x [10.sup.-2]. The reconstruction procedure for a small thin section of the sample is shown schematically in Fig. 2. A 2D image of a section of the full 3D sample is shown in Fig. 1. The interooid vuggy porosities correspond to G(x) = 0 and is automatically introduced in S through Eq. 12 and Eq. 13. Some of the large vugs seen are introduced by the removal of a small fraction of the ooids from the primordial deposit.


Grain decoration

For this model reconstruction, points corresponding to the isopachous layer and the inner cores of the ooids are deposited separately and combined in the following way. First the Nl points in the isopachous layer are decorated with equilateral rhombohedra ([alpha] = -15). The volume of each dolomite [G.sub.i] equals the volume of the associated excluded volume sphere. Although the sphere overlap does not imply crystallite overlap, in this sample the matrix connectivity is retained due to the large overlap and packing . The orientations 6j are chosen randomly from [-20,20] degrees. From the points deposited in inner cores of ooids, we only chose the points [x.sub.i] [??] [M.sub.1] where [M.sub.1] = ([G.sub.1] [union] ([G.sub.2]...... [union]([G.sub.N.sub.1]. The matrix space M of the reconstructed rock sample is fully characterized by a list of N = [N.sub.1] + [N.sub.2] [??] 4 x [10.sup.7] deposited rhombohedral crystallites.

Discretization and property calculation

The reconstructed model is discretized at different resolutions using the 9-point rule described in Eq. 6. Four discretizations are shown in Fig. 3. As in real carbonates, with increasing resolution of discretization, the fraction of the unresolved voxels decreases and porosity increases significantly. Table 1 lists the fraction of the resolved matrix, pore and unresolved voxels for different resolutions of discretization. The rough estimate of the porosity at full resolution obtained by extrapolating these values is in the targeted range 0.25-0.3. The resolved pore scale structure shown in Fig. 3 resembles a typical multiscale microstructure of oolithic carbonate rocks. The strong resolution dependence of permeabilities measured on the discretized samples is also listed in Table 1. Simulation results of resolution dependent pore space and matrix properties of the discretized models have been presented in Biswal et al. (2007). Absolute, single phase permeability and electrical resistivity, or formation factor resistivity, capture the pore space connectivity and geometry, whereas elastic property calculations of bulk and shear moduli recover the connectivity of the matrix structure. These property calculations lend themselves to extrapolation to infinite resolution, like the porosity determined on our model at different discretization levels. This was exemplified in Biswal et al. (2007), and suggests that the investigated properties on multiscale structures can be retrieved consistently by means of scale and resolution dependent simulation from the model.


Here we present the detailed modeling procedure for a second generic carbonate rock microstructure. It is another dolomitized ooid grainstone for which three-dimensional [mu]-CT discretizations are available at three different resolutions. Separate cross-sections for each resolution are shown in Fig. 4. The modeling procedure aims not only to match the morphology but also to incorporate calibrated porosities at different levels of resolution.



Image analysis

Subsample [mu]-CT images of size 1000 x 1000 x 1000 voxels were analyzed for porosity, crystallite types, crystallite size distribution and packing. The greyscale histograms and the threshold dependence of the porosities are shown in Fig. 5 for these three different resolutions. Although the 8.06 [micro]m image seems to have a representative volume for porosity (sidelength of 8.06 mm), the matrix is not well resolved and one observes a more or less continuous histogram. Hence a clear cut-off for the porosity is hard to determine at this resolution. The 1.4 [micro]m data set also seems to have a representative volume (sidelength of 1.4 mm). It appears to resolve the matrix adequately as confirmed by the clear separation of the two peaks in the corresponding greyscale histogram shown in Fig. 5. A similar separation of the two peaks is also observed in the histograms from 0.7 [micro]m resolution images. Ideally, the reference threshold for determining the porosity of the rock should be ascertained from the highest resolution images available. However, in spite of the well resolved matrix, the size of the 0.7 [micro]m images (sidelength of 0.7 [micro]m) is clearly no longer representative. Therefore, to determine the porosity of the rock at different resolutions, we choose a greyscale value of 126 as the reference threshold, determined from the separation of the peaks in the histogram of 1.4 [micro]m image. At this threshold, the 8.06 [micro]m data set has a porosity of 26.27%, the 1.4 [micro]m data set has a porosity of 31.21% and the two data sets with 0.7 [micro]m resolution have porosities of 30.03% and 31.62%, respectively. The large fluctuations observed in the porosities of 0.7 [micro]m images are due to the small volumes of these subsamples. Moreover, the [mu]-CT images at different resolutions analyzed here are cut from different sections of the rock.

In absence of any automated method for estimating the intraooid porosity, rough areal analysis was done for a selection of grains. At the 8.06 [micro]m resolution, no intraooid porosity can be discerned, except in the hollow ooids. Based on visual inspection at different resolutions, approximately 20 percent of the ooids have a hollow inside and around 10-50% of the radius is hollow. The hollow ooids are mostly elongated. The intraooid porosity is found more or less similar at 1.4 [micro]m and 0.7 [micro]m resolutions.


Primordial filter function

A computational model of ooid grain (disc shaped) deposits (Bakke and Oren, 1997) was used as the primordial geometry. In this continuum primordial grain deposits, the grain size distribution was made to match with that extracted from the [mu]-CT samples. The ellipsoidal ooid deposits (Fig. 6, left) are all perfectly oriented. The side length of the cubic box with primordial deposits is 1.5 mm. In order to avoid unwanted border effect it is taken from the center of the deposited cubic box of sidelength 1.8 mm with 2152 ooid (disc shaped) deposits. After gridding at 300 x 300 x 300 voxels, the porosity of the disc grain pack is roughly 34%. This estimate is somewhat lower than the actual porosity of the continuum grainpack. To compensate the intraooid porosity that will be introduced later and still be able to achieve the targeted porosity, the sizes of the primordial grains are enlarged by 2%. The greyscale primordial filter function G(x) is then constructed using Eq. 10.

Based on the image analysis, five different kinds of grain parameters (primordial phases) were chosen. The primordial phase 1 corresponds to the isopachous rim on all the grains. The remaining four primordial phases are defined as follows. Let [H.sup.j.sub.i] denote the ith primordial grain in the jth primordial phase ({[H.sup.j.sub.i], i = 1,2,..,[m.sub.j], j = 2,..,5}). [m.sub.j] is the number of grains in primordial phase j and m = [summation.sup.5.sub.i=2] [m.sub.j] is the total number of primordial grains. 30% of the primordial grains are randomly selected and defined as primordial phase 2. They are denoted by {[H.sup.2.sub.1], i = 1, 2,..,[m.sup.2]} and correspond to the tight intraooid crystallite packing. 20% of the randomly chosen primordial grains are grouped into the primordial phase 3 corresponding to the moderate intraooid crystallite packing and 30% of the randomly chosen grains are grouped into the primordial phase 4 corresponding to the loose intraooid crystallite packing. The remaining 20% of the primordial grains are grouped into the primordial phase 5 that contains grains with tight intraooid crystallite packing along with hollow, vuggy cores. We define new primordial filter functions separately for each of the primordial phases as


where G(x) is defined in Eq. 10.

Deposition of centers and grain parameters

Model parameters for different primordial phases are chosen as follows and points were deposited randomly and sequentially in S. For computational efficiency we divide S [subset] [R.sup.3] into smaller nonoverlapping cubic cells U each of sidelength l, i. e., S = [U.sub.1][union][U.sub.2][union][U.sub.3]...Points are deposited randomly in each cell [U.sub.i] and added to the final list if the overlap condition defined in Eq. 3 is satisfied.

1. For primordial phase 1 (isopachous rims of all the ooids with moderate packing of intraooid crystallites), we chose [R.sub.min] = 5 [micro]m, [R.sub.max] = 12.5 [micro]m, [lambda] = 0.6, l = 37.5 [micro]m and [rho] = 1.8335 x [10.sup.-3] and set Pr([[OMEGA].sub.1]) = 0, for

[[OMEGA].sub.1] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, [H.sup.1]([x.sub.i]) [??] (0,0.15)} (15)

2. For primordial phase 2 (30% of the randomly chosen ooids with tightly packed intraooid crystallites) we chose [R.sub.min,] = 3.75 [micro]m, [R.sub.max] = 11.25 [micro]m, [lambda] = 0.8, l = 33.75 [micro]m and [rho] = 5.0301 x [10.sup.-3] and set Pr([[OMEGA].sub.2]) = 0, for

[[OMEGA].sub.2] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, [H.sup.2] ([x.sub.i])[member of][0.0.15]}. (16)

3. For primordial phase 3 (20% of the randomly chosen ooids with intermediate packing, i. e., less than the moderate packing in the rim,) we chose [R.sub.min] = 3.75 [micro]m, [R.sub.max] = 11.25 [micro]m, [lambda] = 0.5, l = 33.75 [micro]m and [rho] = 2.012 x [10.sup.-3] and set Pr([[OMEGA].sub.3]) = 0, for

[[OMEGA].sub.3] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, [H.sup.3] ([x.sub.i])[member of][0.0.15]}. (17)

4. For primordial phase 4 (30% of the randomly chosen ooids with loosely packed intraooid crystallites) we chose [R.sub.min] = 3.75 [micro], [R.sub.max] = 11.25 [micro]m, [lambda] = 0.5, l = 33.75 [micro]m and [rho] = 1.6767 x [10.sup.-3] and set Pr([[OMEGA].sub.4]) = 0, for

[[OMEGA].sub.4] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, [H.sup.4] ([x.sub.i])[member of][0.0.15]}. (18)

5. Intraooid vugs are introduced in the remaining 20% of the ooids deposited with parameters [R.sub.min] = 3.75 [micro]m, [R.sub.max] = 11.25 [micro]m, [lambda] = 0.8, l = 33.75 [micro]m and [rho] = 5.0301 x [10.sup.-3]. Vug radius varies randomly (with uniform distribution) from 10% to 50% of the corresponding ooid radius. With each of the grains in the primordial phase 5, we associate a random number {[[eta].sub.i], i [less than or equal to] [m.sub.5]} that has a fixed but random value in [0.5,1]. We set Pr([[OMEGA].sub.5]) = 0, for

[[OMEGA].sub.5] = {[omega][member of][OMEGA]:[there exists]0[less than or equal to]i[less than or equal to]N, [H.sup.5] ([x.sub.i])[??][[eta]([x.sub.i]),0.85]}, (19)

where [[eta].sub.5]([x.sub.i]) = [[eta].sub.j] if [x.sub.i] [member of] [H.sup.5.sub.j] for j [less than or equal to] [m.sub.5]

The points are then decorated with equilateral rhombohedra. Each rhombohedron [G.sub.i] (before rotation) centered at [x.sub.i] is the intersection of three pairs of parallel planes as defined in Eq. 7. The distance [d.sub.i] is chosen such that the volume of [G.sub.i] equals the volume of the associated excluded volume sphere and [alpha] = -15 degrees. Each rotation angle [[theta].sub.j] is chosen randomly from [-60,60] degrees. Roughly 6.4 x [10.sup.6] rhombohedral crystallites are sufficient to fill the whole 1.5 mm x 1.5 mm x 1.5 mm cubic sample. A cross-section of the discretized sample with these decorated crystallites is shown in Fig. 6.



The reconstructed samples are then discretized with 125 collocation points (corresponding to n = 5) as described above in the section entitled "SYNTHETIC [mu]-CT". Seven discretized data sets are generated. One data set contains the full sample discretized at 8.06 [micro]m. Next, two non-overlapping subsamples, half the size of the full sample, were discretized at 1.4 [micro]m resolution. These two subsamples are the lower left corner and the upper right backwards corner of the full sample, i. e., the second subsample is displaced along the diagonal. Finally, there are four non-overlapping subsamples cut along the diagonal, each one quarter the size of the full sample and discretized at 0.7 [micro]m resolution. Cross-sections of these greyscale discretizations are visualized in Fig. 4 for all three resolutions. Porosity was measured at different thresholds for all these discretized samples. Threshold dependent porosities averaged over the subsamples are presented in Table 2. Porosities of the [mu]-CT images thresholded at 126 (see section "SYNTHETIC [mu]-CT") were compared with model discretizations at corresponding resolutions for threshold values of [m.sub.c] = 48, 52, and 56. The porosities of the model discretizations correspond well to those of the [mu]-CT images at [m.sub.c] = 52. Like the 0.7 [micro]m resolution ,u-CT images large porosity fluctuations are also observed in the discretized samples at this resolution. Morphologically, the discretized samples shown in Fig. 4 also match well with the [mu]-CT images shown in Fig. 4.

It is important to note that the above discretizations represent synthetic [mu]-CT images from the model and owing to the continuum representation, can be generated at arbitrary resolutions. This allows investigation of the rock properties at both intermediate and higher resolutions than the available [mu]-CT images. However, a direct comparability with experimental [mu]-CT images requires addition of a small amount of Gaussian white noise to these synthetic images. For example, to match the greyscale densities, we first shift the peaks and width of the greyscale density histograms from the model to match with that of the the experimental [mu]-CT images. Then we add a small amount of noise (see Fig. 7) to match the peak corresponding to the matrix phase. The discretizations at resolutions 8.06 [micro]m, 1.4 [micro]m and 0.7 [micro]m are added with white noise of standard deviation 19, 17 and 15 respectively. These images are shown in Fig. 7 and the threshold dependent porosity reveals an impressive match with that measured from the experimental [mu]-CT images as shown in Fig. 5.

Two-point correlation function

Quantitative characterization and comparison of the microstructure is carried out through the two-point correlation function (Oren et al., 2007) measured from the discretizations of the model and the original [mu]-CT images at different resolutions. These greyscale representations are thresholded to obtain a digitized representation I(x). The indicator function I(x) = 1 if x [member of] matrix space and I(x) = 0 if x [member of] pore space, where x is the position vector of the voxel on the digitized grid. Assuming homogeneity, the two-point correlation function C(r) is defined as

C(r) = <(I(x) - [psi])(I(x+r) - [psi])>/<[(I(x) - [psi]).sup.2]), (20)

where [psi] is the porosity of the sample. For the following analysis C(r) is computed for r aligned in the X, Y and Z directions.


For this analysis, the [mu]-CT images ([1000.sup.3] voxels, 8-bit greyscale) at resolutions 8.06 [micro]m, 1.4 [micro]m and 0.7 [micro]m are digitized at threshold value of 126. The model discretizations at resolutions 8.06 l.un (full sample, [186.sup.3] voxels), 1.4 [micro]m (two subsamples, 5363 voxels) and 0.7 [micro]m (four subsamples, [536.sup.3] voxels) are digitized at threshold value [m.sub.c] = 52. The two-point correlation functions measured on these digitized representations are plotted in Fig. 8. At 8.06 [micro]m resolution, the model measurements agree reasonably well with the [mu]-CT data in the short range with marginal discrepancy in the long range. At 1.4 [micro]m resolution, the model measurements also agree well with the [mu]-CT data. At 0.7 [micro]m resolution, model measurements from several of the subsamples agree with the [mu]-CT data in the short range. The nonrepresentative sample sizes give rise to large sample to sample fluctuations making a comparison at this resolution difficult. The slope of the tangent drawn at C(0) provides a rough estimate of the specific surface area of the solid/void and is related to the permeability of the sample (Hilfer, 1996). The above results indicate that along with the quantitative agreement in the porosities (Table 2) and in the morphologies (Fig. 4 and Fig. 7), the transport properties measured from the model subsamples at different resolutions are also expected to show reasonable agreement with the subsamples of the rock used for the [mu]-CT imaging.

Discrepancies observed between the [mu]-CT images and the model discretizations may originate from a number of different sources. The model parameters are only approximate. The dolomite sizes and their overlap are measured from two dimensional micrographs. The phases are determined from rough areal analysis of some selected grains in the experimental [mu]-CT images shown in Fig. 4. The primordial ooid deposits differ from the ooids visible in the [mu]-CT images in a number of ways. The primordial ooids have circular cross-sections in the Z-direction and elliptical cross-section in the other two directions, thus the shape of the primordial ooids match with the original rock only in the X and Y directions. As a result of the deposition procedure (Bakke and Oren, 1997), they are perfectly aligned with respect to Z-direction, whereas the orientation of the ooids in the original rock is more isotropic. Due to these reasons a large discrepancy is observed between the correlation functions measured in the Z-direction and the other two directions for the model (Fig. 8). This emphasizes the need for using valid primordial input and model parameters in the model reconstruction to achieve good correspondence with the spatial correlations in the microstructure and transport parameters.


In this paper we present a model for multiscale porous media that provides a continuum description of the porous rock at arbitrary precision. It reproduces the enormous variation of "pore sizes" (vugs, pores, cracks, slits etc) that is typical of multiscale porous media. The range of length scales that can be included is limited only by the floating point precision of computers. The petrophysical properties of the reconstructed models can be studied as a function of resolution without any limit in the magnification, i.e., scale dependence of porosity, connectivity or network models. The model permits to represent vugginess (both separate, moldic and touching such as Breccia in the sense of Lucia (1999)). It reproduces intergranular porosity, intragranular porosity, correlation with primordial depositional texture and full connectivity of pore and matrix space. The model parameters allow calibration of porosities and physical parameters at different resolutions.

A successful reconstruction of a given rock using this model depends on the accuracy of both the input parameters obtained from image analysis and the primordial filter function. The matrix connectivity is ensured by finite crystallite overlap that is implemented through an overlap of spheres associated with each crystallite. Although it makes the algorithm simple and fast, this approach may not be suitable for rocks containing elongated or complex objects. To ensure faster and computationally feasible deposition of millions of spheres, the algorithm presented in this paper uses a subdivision of the rock into sub-cells, filled sequentially. A faster algorithm for depositing hundreds of millions of objects with a given size distribution, satisfying simultaneously the overlap rule, would be necessary to overcome these restrictions. To the best of our knowledge, the proposed modeling approach is currently the only feasible method available for generating complex multiscale carbonate pore-scale models of laboratory scale rock samples.

This process inspired continuum model represents a feasible pore scale modeling technology for multiscale porous media such as carbonates. Many elements of realistic diagenetic processes are included in the modeling procedure. The model was tested on two examples of oolithic dolostones and can be easily generalized to model a wide variety of carbonate rocks. It can also be used to reconstruct other kinds of multiscale porous media such as sandstones with strong heterogeneities and fractured porous media.


BB and RH gratefully acknowledge funding of this work by StatoilHydro ASA, Norway.

(Accepted January 30, 2009)


Adler PM (1992). Porous Media--Geometry and Transports. Boston: Butterworth-Heinemann.

Anselmetti FS, Luthi SM, Eberli GP (1998). Quantitative characterization of carbonate pore systems by digital image analysis. AAPG Bulletin 82:1815-36.

Arns CH, Knackstedt MA, Mecke KR (2003). Reconstructing complex materials via effective grain shapes. Phys Rev Lett 91:215506.

Bakke S, Oren PE (1997). 3-d pore-scale modeling of sandstones and flow simulations in pore networks. SPE J 2:136-49.

Biswal B, Oren PE, Bakke S, Held RJ, Hilfer R (2007). Stochastic multiscale model for carbonate rocks. Phys Rev E 75:061303.

Dunham R (1962). Classification of carbonate rocks according to depositional texture. In: W. Ham, eds. Classification of carbonate rocks: American Association of Petroleum Geologists, 108-21.

Feder J (1980). Random sequential adsorption. J Theor Biol 87:237-54.

Fernandes CP, Magnani FS, Philippi PC, Ddian JF (1996). Multiscale geometrical reconstruction of porous structures. Phys Rev E 54:1734-41.

Hilfer R (1996). Transport and relaxation phenomena in porous media. Adv Chem Phys 92:299-424.

Hilfer R (2000). Local porosity theory and stochastic reconstruction for porous media. In: D. Stoyan, K. Mecke, eds. Rdumliche Statistik and Statistische Physik. Lect Notes Phys 254. Berlin: Springer, 203-41.

Lucia FJ (1999). Carbonate Reservoir Characterization. Berlin: Springer.

Moore CH (1989). Carbonate Diagenesis and Porosity. Amsterdam: Elsevier.

Moore CH (2001). Carboante Reservoirs - Porosity Evolution and Diagenesis in a Sequence Stratigraphic Framework. New York: Elsevier.

Okabe H, Blunt MJ (2007). Pore space reconstruction of vuggy carbonates using microtomography and multiple-point statistics. Water Resour Res 43:W 12S02.

Oren PE, Bakke S, Held RJ (2007). Direct pore-scale computation of material and transport properties for north sea reservoir rocks. Water Resour Res 43: W 12504.

Song YQ, Rye S, Sen P (2000). Determining multiple length scales in rocks. Nature 406:178-81.

Stauffer D, Aharony A (1992). Introduction to Percolation Theory. London: Taylor and Francis.

Stoyan D, Kendall WS, Mecke J (1995). Stochastic Geometry and Its Applications. Chichester: John Wiley & Sons.

Thovert JF, Yousefian F, Spanne P, Jacquin CG, Adler PM (2001). Grain reconstruction of porous media: Application to a low-porosity fontainebleau sandstone. Phys Rev E 63:061307.


(1) ICP, Universitat Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart, Germany, (2) S. V College, University of Delhi, New Delhi--110 021, India, (3) Numerical Rocks AS, N-7041 Trondheim, Norway, (4) StatoilHydro ASA, N-7005 Trondheim, Norway, (5) Institut fur Physik, Universitdt Mainz, 55099 Mainz, Germany e-mail:,,,,
Table 1. Fraction of resolved pore voxels ([f.sub.p]), unresolved
voxels ([f.sub.u]) and matrix voxels ([f.sub.m]) at different levels of
discretization. In column 5 and 6 the absolute permeabilities computed
on the discretized samples are listed. To compute the permeability, the
unresolved voxels are either converted to matrix ([k.sub.m]) or to pore

(in [micro]m) [f.sub.p] [f.sub.u] [f.sub.m]

20. 0.12683 0.8463 0.02687
13.333 0.13868 0.78332 0.078
10. 0.14595 0.7154 0.13864
6.6667 0.15562 0.58397 0.2604
5. 0.1628 0.4756 0.3616
4. 0.16914 0.3941 0.43677
3.333 0.1749 0.3341 0.491
2.5 0.18506 0.254 0.561
0 (extrapolated) 0.253 0.008 0.739

(in [micro]m) [K.sub.m] (in mD)[K.sub.p]
(in mD)

20. 256.32 1.8467 x [10.sup.7]
13.333 1643.2 3.7239 x [10.sup.6]
10. 3495. 1.0391 x [10.sup.6]
6.6667 6539.5 2.0192 x [10.sup.5]
5. 8150.4 7.3075 x [10.sup.4]
4. 8520.3 4.0128 x [10.sup.4]
3.333 8570.9 2.7383 x [10.sup.4]
2.5 8558.6 1.834 x [10.sup.4]
0 (extrapolated)

Table 2. Variation in threshold dependent porosities of the [mu]-CT
images and the discretized samples (without noise) at three different
resolutions. [mu]-CT measurement for 0.7 [micro]m resolution is
averaged over two different subsamples. Model measurements for 1.4
[micro]m and 0.7 [micro]m resolutions are averaged over two and four
subsamples respectively.

Resolution Original [phi] (n) Model
(in [micro]m) [mu]-CT n = 126

8.06 mCT_F 0.2627 model_F
1.40 mCT_H 0.3121 model_l
0.70 mCT_Q 0.3082 model_Q

Resolution [phi] ([m.sub.c])
(in [micro]m)
 [m.sub.c] = 48 [m.sub.c] = 52 [m.sub.c] = 56

8.06 0.2537 0.2615 0.27
1.40 0.3096 0.3133 0.317
0.70 0.3077 0.3097 0.3116
COPYRIGHT 2009 Slovenian Society for Stereology and Quantitative Image Analysis
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2009 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Biswal, Bibhu; Oren, Pal-Eric; Held, Rudolf J.; Bakke, Stig; Hilfer, Rudolf
Publication:Image Analysis and Stereology
Article Type:Report
Geographic Code:4EUGE
Date:Mar 1, 2009
Previous Article:A stereological unfolding method for the study of the mitochondrial network.
Next Article:Velocity field computation in vibrated granular media using an optical flow based multiscale image analysis method.

Related Articles
Model optimizes moisture transport.
Modeling and numerical analysis of the two-phase fluid flow through porous medium.
Diffusion-thermo and thermal-diffusion effects on free convective heat and mass transfer flow in a porous medium with time dependent temperature and...
Carbon nanotube devices; properties, modeling, integration and applications.
Porous media; heat and mass transfer, transport and mechanics.
Nano- and micromechanics of polymer blends and composites.
Advances in strength of materials.
3D morphological modelling of a random fibrous network.
Effect of thermal diffusion on an oscillatory MHD convective flow with mass transfer through a porous medium.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters