# Fractal geometry of sedimentary rocks: simulation in 3-D using a Relaxed Bidisperse Ballistic Deposition Model

Abhra Giri, Sujata Tarafdar, Philippe Gouze, Tapati Dutta

► To cite this version:

Abhra Giri, Sujata Tarafdar, Philippe Gouze, Tapati Dutta. Fractal geometry of sedimentary rocks: simulation in 3-D using a Relaxed Bidisperse Ballistic Deposition Model. Geophysical Journal International, 2013, 192 (3), pp.1059-1069. 10.1093/gji/ggs084. hal-00823061

HAL Id: hal-00823061

https://hal.science/hal-00823061

Submitted on 11 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Geophys. J. Int. (2013) 192, 1059–1069

doi: 10.1093/gji/ggs084

# Fractal geometry of sedimentary rocks: simulation in 3-D using a Relaxed Bidisperse Ballistic Deposition Model

Abhra Giri,<fnref n="1" />,<fnref n="2" /> Sujata Tarafdar,<fnref n="2" /> Philippe Gouze<fnref n="3" /> and Tapati Dutta<fnref n="1" />

Accepted 2012 November 22. Received 2012 November 21; in original form 2012 June 17

## SUMMARY

Several studies, both theoretical and experimental, show that sedimentary rocks have a fractal pore-grain interface. In this paper a computer simulated 3-D sedimentary rock structure generated by the Relaxed Ballistic Bidisperse Deposition Model (RBBDM), is investigated to characterize the micro structure of its pores. The pore volume and the rock-pore interface show the same fractal dimension indicating that the pore volume is a fractal. The two point density correlation is computed for the pore space and the results compare favourably with the range reported from experiments. An array of 2-D X-ray tomography micrograph sections of a real sedimentary rock, an oolitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France), was used to generate a 3-D bitmap. The 3-D real rock sample generated in this manner, was analysed for similar studies as the simulated structure. The results were compared with those obtained from simulation. The simulation results agree qualitatively with the real rock sample. Diffusion through the connected pore space of the simulated structure was studied using a random walk algorithm and the results compared with the similar simulation study done on the 3-D oolitic limestone specimen. In both cases diffusion was found to be anomalous indicating that the sedimentary rock has a fractal geometry. The favourable comparability of results between the simulated and real rock supports the usefulness of the model of sedimentary rock generation which can be applicable to transport phenomena.

**Key words:** Image processing; Fractals and multifractals; Microstructures; Sedimentary basin processes; mechanics, theory and modelling.

## 1 INTRODUCTION

Transport of fluids through the connected pore space in sedimentary rocks is of profound importance in the extraction process of oil and natural gas, the problem of CO₂ sequestration, groundwater flow, to give a few examples. This transport process is intimately connected to the nature of the pore space which often gets modified if the flow is reactive diffusive and/or advective. In this work we intend to investigate and characterize the pore and rock phase of the sedimentary rock generated by the RBBDM. Several studies on sedimentary rocks (Katz 1985; Wong 1987; Krohn 1988a,b) have suggested that the pore system has a fractal nature. Although the pore-rock interface has been shown definitely to be a fractal, there is a controversy about whether the pore space is a volume fractal too. Not many measurements have been possible to probe the pore volume directly. Any fractal nature of the pore space is possibly a consequence of the formation of the rocks by sedimentation and diagenesis of granular matter.

In an earlier work (Roy & Tarafdar 1998; Dutta & Tarafdar 2003), the authors simulated a 3-D sedimentary rock structure using the

Bidisperse Ballistic Deposition Model (BBDM) and investigated its pore structure. The model mimicked the deposition process of grains under gravity to generate the rock structure. This model was extended under the 'Relaxed Bidisperse Ballistic Deposition Model' (RBBDM; Sadhukhan et al. 2007) where the deposited grains were allowed to relax to mimic the process of compaction. Although the RBBDM has been used to study transport properties such as permeability and conductivity through the pore space (Sadhukhan et al. 2007a,b, 2008), the pore structure generated by this process has not been investigated. The study of the microstructure is very important if any model for sedimentary rock generation is to be used for the investigation of transport properties. We shall investigate the effect of the relaxation of the grains that mimics the process of compaction, on the micro geometry of the pore space by comparing our results with that obtained in (Dutta & Tarafdar 2003).

In the following sections we briefly discuss the model RBBDM used for generating the sedimentary rock structure of varying porosity. Two different grain sizes are used in the deposition process. The dependence of porosity on the fractional composition of the grains is studied. The effect of relaxation of the grains is studied with

Downloaded from https://academic.oup.com/gji/article/192/3/1059/819491 by guest on 11 June 2021

<footnote n="1">¹Physics Department, St. Xavier's College, Kolkata 700016, India. E-mail: tapati_mithu@yahoo.com</footnote>

<footnote n="2">²Condensed Matter Physics Research Centre, Physics Department, Jadavpur University, Kolkata 700032, India</footnote>

<footnote n="3">³Geosciences, Université de Montpellier 2, CNRS, Montpellier, France</footnote>
respect to the BBDM model. Two-point density correlation function is calculated for the entire pore space of the generated rock structure. The results are compared with the values calculated on real rock samples and those obtained from BBDM (Dutta & Tarafdar 2003). A study of the nature of both, the connected pore space and the total pore space in the generated rock structure, is done. In the study using BBDM (Dutta & Tarafdar 2003), we had observed that only the connected pore space has a fractal nature. If the total pore space which is made up of all isolated pores and also the connected pore is considered, there is no evidence of fractal behaviour. The contribution from the disconnected pores masks any signature of a fractal nature. The pore space becomes homogeneous on average. In this study we find that the behaviour in both cases is exactly the same, that is, they exhibit the same fractal nature. The rock-pore interface is also checked for any fractal nature. The results are compared with those of BBDM to study the role of compaction of grains. These results are also compared with the calculations done on real sedimentary rocks. We also study diffusion of a random walker in the connected pore space. A random walk through a fractal pore volume is expected to show anomalous behaviour. This is evident from our study and reinforces the fractal picture of the pore space. Finally, we draw the conclusions of our investigation and share some of the future plans on this subject.

## 2 MODEL

The details of the BBDM and RBBDM are discussed in Dutta & Tarafdar (2003) and Sadhukhan *et al.* (2007), respectively. However the basics of the RBBDM will be described briefly here for the sake of completeness. The porous structure is generated by ballistic deposition of grains of two different sizes. We drop cubes 1 × 1 × 1 and parallelepiped 2 × 1 × 1 grains on a square substrate. Models on disordered systems usually consider particles to be spherical, unlike RBBDM. However real sand grains are often angular (Pettijohn 1984). In this respect our approximation of rectangular parallelepiped grains is no worse than a spherical grain approximation. Due to abrasion at points of contact, mature grains meet along planes and tend to show a slight elongation. A ratio of long to short axis of 1.0–2.5, generally close to 1.5, is not uncommon. So the aspect ratio 2 is realistic. The cubes are chosen with a probability *p* and elongated grains with probability (1 − *p*). The presence of the longer grains leads to gaps in the structure. The porosity *φ*, defined as the vacant fraction of the total volume, depends on the value of *p*.

For *p* = 1, a compact rock structure with zero porosity is formed. As *p* is decreased, isolated pore ‘clusters’ start appearing and the porosity increases. For a specific value of *p*, a structure spanning cluster is generated. However this cannot be called a ‘percolation threshold’ as in the case of random percolation (Stauffer & Aharony 1994). The RBBDM, being a modification of the Random Deposition Model, for *p* = 1, the surface width keeps increasing with height. In fact, in the limit of infinite height, at least one narrow structure spanning pore cluster is always present. Thus one cannot speak of a ‘percolation threshold’. When *p* is gradually decreased to below 1, larger grains are introduced and the grains settle on the structure following the Ballistic Deposition Model. The introduction of even an infinitesimal quantity of the larger grain, can close a deep surface trench creating an elongated pore cluster. Obviously these pore clusters are longer near the surface of the structure than at the bottom. The presence of a single large grain sitting atop a long pore cluster, introduces correlation between adjacent columns. As

the fraction of large grains increase, the correlation spreads through the system. So a substrate of sufficient height needs to be generated so that the porosity value can stabilize. This has been demonstrated for the BBDM (Karmakar *et al.* 2005). In fact in spite of relaxation of longer grains in RBBDM, the situation here is similar to BBDM. Our model is different from the random percolation problem as even an infinitesimally small fraction of larger grains introduces correlation between columns, thus robbing the system of its randomness.

The RBBDM has the potential of generating a structure with a connected rock phase that is needed for any stable structure, and a tunable porosity. As the fraction of longer grains is increased, unstable overhangs can develop. If a larger particle settles on a smaller particle, a one-step overhang is created. If a second larger particle settles midway on the previous large particle, a two-step overhang is created, when there is no supporting particle immediately below the protrusion of the second overhang. This two-step overhang is not stable and the second large particle topples over if possible, according to the rule shown schematically in Fig. 1. This leads to compaction.

It had been discussed in BBDM (Manna *et al.* 2002; Dutta & Tarafdar 2003) that the sample attains a constant porosity only after a sufficient number of grains (depending on sample size) have been deposited to overcome substrate effects. Here, a $Lx \times Lx \times Lz$ size sample was generated, from which a $Lx \times Lx \times Lx$ sample was selected after the porosity had stabilized to within a fluctuation of 0.001 per cent. The selected sample was chosen from below the deepest trough at the surface to eliminate surface effects. All simulations were carried out on this sample. To check for finite size effects, we carried out our studies for $Lx = 32, 64, 128, 256$ for which $Lz = 1000, 2000, 4000, 7000$, respectively. The results reported in this work did not show any finite size dependence. All results on simulation are reported for $256 \times 256 \times 256$. An elongated grain is expected to have its long axis vertical, while falling under gravity. However, when it reaches the ground, its stable position is more likely to be with the long axis horizontal as this will keep the centre of mass of the grain at a lower height. This will in turn lower the total energy and stabilize the structure. This process was mimicked in the deposition algorithm where the elongated grains are deposited with their long axis horizontal and parallel to either the x- or y-axes with equal probability. The vertical direction is the direction of the z-axis and coincides with the direction of grain deposition.

Sedimentary rocks usually have a porosity that rarely exceeds 0.5. As the fraction of larger grains is increased, the porosity of the sample increases. A vertical (x-z), and horizontal section (x-y), of the sample at its maximum porosity $\phi_{max} = 0.45$ are shown in Figs 2(a) and (b), respectively. It is evident that the sample looks

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Toppling rule of the larger grains—when a larger grain develops a two-step overhang, marked 1 in the figure, with at least two vacant sites immediately below the overhang, it topples over in the direction indicated by the arrow to assume a more stable state, marked 2.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Panel (a) shows *x–z* section of simulated structure for *φ*<sub>max</sub> = 0.45, that is, high porosity. (b) *x–y* section at same porosity. Structure looks more isotropic. The white indicate pore clusters.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Panels (a) and (b) show *x–z* and *y–z* section of simulated structure for *φ* = 0.072, that is, low porosity, matching the porosity of the real rock. The *z*-axis indicates the vertical direction. The white indicate pore clusters. Panel (c) *x–y* section at same porosity. Anisotropy in pore cluster structure is quite pronounced. The white indicate pore clusters.</figcaption></figure>

quite isotropic and starts to resemble a real rock. Vertical sections (*x–z* plane) and (*y–z*) planes, of the generated sample at very low porosity *φ* = 0.072, are shown in Figs 3(a) and (b). A horizontal section (*x–y* plane) of the sample at the same porosity value is shown in Fig. 3(c). The anisotropy in the pore geometry is clearly visible. At low porosities, elongated and isolated pore clusters appear along the *z*-direction, the direction of particle assembly, while the distribution of pores along the horizontal plane is quite homogeneous.

In our model, the presence of large grains introduces correlation between neighbouring grains. The toppling rule of the larger grains decreases the porosity somewhat. Therefore as the fraction of large grains increase, the porosity increases though not monotonically as in the BBDM. Compaction of the grains results in the porosity showing a maximum at *p* = 0.5. To compare our simulation results with real rock samples, X-ray tomography micrographs of 2-D sections of real sedimentary rock sample obtained from an oolitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). This rock has been used by (Luquot &

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Each section is a square of side 2.58 × 10⁻³ m. Panels (a) and (b) show sections of the real rock cut in mutually perpendicular planes. These are perpendicular to bedding planes. Panel (c) shows a section of the bedding plane. These are sections of oolitic limestone(pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). The graphical files were constructed from X-ray microtomographs of 2-D sections of the rock. The real structure seems to be isotropic.</figcaption></figure>

Gouze 2009) and is essentially composed of recrystallized ooliths with a mean diameter of less than a few hundred micrometres. Each pixel of the micrograph corresponds to 5.06 µm. Each section was converted to a binary file form such that 0 corresponded to a pore site and 1 corresponded to a rock site. A binary file picture of the real rock is obtained by suitable thresholding of a grey-scaled picture, using MATLAB, Fig. 4. An array of 1000 consecutive sections were put together precisely to reconstruct the binary file form of the real 3-D rock structure. The real structure was 1000 × 1000 × 1000 voxels in size, and all study on real rock was carried out on this structure. 2-D sections, (*x–z*) plane and (*y–z*) plane, cut in the direction of assembly (growth) of the reconstructed 3-D structure, are shown in Figs 4(a) and (b), respectively. Comparison between Figs 4(a) and (b) shows the pore distribution to be isotropic. A 2-D section of the bedding plane of the rock structure is shown in Fig. 4(c). There seems to be slight anisotropy in the pore distribution in the bedding plane and the direction of growth.

For every *p* value studied, we have calculated the total porosity of the sample and also the porosity of the connected pore space. The variation of porosity with the fraction of cubes *p* obtained in the RBBDM, for both these cases is shown in Fig. 5. The two graphs are almost parallel to each other. The porosity shows a maximum at *p* = 0.5 unlike BBDM (Dutta & Tarafdar 2003), where a monotonic decrease of porosity with increase in *p* was noted. From Fig. 5 it appears that there are some *p* values for which combinations of larger grains to smaller grains yield the same value of porosity. A case study of *φ* = 0.44 may be made in this respect. It is well known that sedimentary rocks that have same value of porosity, may have very different pore structure leading to very different transport properties. The shape and size distribution of pores and their connectedness is responsible for such variations. In our model, the same value of porosity *φ* = 0.44, for example, is obtained for *p* = 0.2 and also *p* = 0.7. The difference lies in the micro geometry of the pore space. This point is dealt with in detail in the next section where we have studied the two-point density correlation function
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Variation of porosity with *p*. The solid line denotes the total porosity of the simulated structure, while the dotted line corresponds to the porosity of the connected pore space. The values corresponding to the real rock studied are marked on each graph.</figcaption></figure>

that throws light on the pore distribution. The RBBDM produces a porosity range from $\phi_{max} = 0.451-0.0$. The increase in fraction of large grains has two competing effects: (1) the creation of more gaps increasing porosity and (2) the increased toppling events tending to decrease $\phi$. This competition leads to a non-monotonic variation in $\phi$ with $p$ with a maximum porosity $\phi_{max} = 0.45$.

The total porosity of the real rock studied was 0.073 which corresponds to a $p = 0.995$, Fig. 5. The porosity of the connected pore space of the real sample was found to be 0.067 which corresponds to $p = 0.970$. From Fig. 5 we see that a connected pore system is obtained upto $p = 0.99$ for the finite system we have considered in the simulation. The number of pore sites (voxels) at this value of $p$ is 1532 which is greater than $Lx$, the system size. Clearly the connected pore cluster at this stage is not a straight channel but has a tortuous path.

The RBBDM can create any porosity in the range $\phi_{max}$ to 0 by varying $p$. This covers the range of porosity of real sedimentary rocks which have a porosity that rarely exceeds 0.5. The real sample we have analysed has $\phi = 0.073$ ($p = 0.997$) which falls in the range of RBBDM, but the model gives more realistic structures for higher $\phi$. Very low, $\phi \sim 0.1$ or less, creates mostly long linear pores as seen in Fig. 3.

## 3 TWO-POINT DENSITY CORRELATION FUNCTION

The two-point density correlation function is defined as the expectation value that two points are separated by $\mathbf{r}$ and belong to the same phase, though not necessarily connected by the same phase. Mathematically,

$$S_2 = \langle \rho(\mathbf{r}')\rho(\mathbf{r} + \mathbf{r}) \rangle, \quad (1)$$

where the average is done over the entire pore volume and all values of $\mathbf{r}'$ are considered. $\rho$ takes value 1 if there exists a pore, else it takes the value 0. Ordinary fractals are typically isotropic implying that the density correlation depends only on $\mathbf{r}$ and not on the direction.

In our calculation, the two-point density correlation has been calculated for the pore phase in all the three orthogonal directions

and for every pair of points $\mathbf{r}'$ and $(\mathbf{r}' + \mathbf{r}_e)$ with $r = 0, 1, 2, ..., r_c$. $r_c$ is the cut-off value determined by system size, $\mathbf{e}_i$ is the unit vector along the $i$th orthogonal axis. Therefore a calculation of the two-point density correlation function along the $z$-axis involves

$$S_2(z) = \langle \rho(z')\rho(z' + z) \rangle. \quad (2)$$

A study of $S_2(x)$, $S_2(y)$ and $S_2(z)$ reveals transverse isotropy in the $x-y$ plane with $S_2(x) = S_2(y)$. Their value is different from $S_2(z)$ showing anisotropy along this direction which is attributed to the structure growth process.

The results of our simulation resemble the nature of Berea sandstone, Ironton-Gallesville sandstone (Berryman & Blair 1986) and the Fontainbleau sandstones as reported by (Manwart et al. 2000), none of which have significant clay content. $S_2(x)$ and $S_2(z)$ have been plotted against $\mathbf{r}$ at the maximum porosity $\phi = 0.457$, and a low value of porosity $\phi = 0.352$ in Fig. 6(a). Thus the impenetrable grains of our simulation seem to have a similarity to clean sandstone samples. The dips in the curves give information of the average grain cluster size. A larger dip for low porosity values compared to high porosity along the $z$-direction is evident from Fig. 6(a). The inset in Fig. 6(a) shows the same study at very low porosity $\phi = 0.077$ which is almost the same as the real rock sample studied.

The variation of the two-point density correlation function for $p = 0.2$ and $0.7$, both of which have the same value of porosity $\phi = 0.44$, is depicted in Fig. 6(b). For lower $p$ value, that is, 0.2, the minimum is higher and occurs at slightly shorter values of $\mathbf{r}$ than at the higher $p = 0.7$. This is highlighted in the inset of Fig. 6(b). The lower minimum of the two-point density correlation function for $p = 0.7$ indicates greater anti-correlation between pores implying that the distribution is more ordered here. The shorter distance of the minimum at $p = 0.2$ implies that the average grain cluster is smaller in size than at $p = 0.7$. These differences though small, are unmistakable. Thus the non-unique combinations of $p$ that exist for some values of porosity is due to variation in the shape and size of the pore clusters. The anisotropy in pore distribution is clearly evident at lower values of porosity and the structure becomes more isotropic at higher porosities. The minimum in the density correlation function occurs at greater distance for low porosity and shifts
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Two-point density correlation along x- and z-direction of simulated rock structure for high and low porosity values. The anisotropy decreases as porosity increases. The inset shows the variation for $\phi = 0.077$. Panel (b) shows the variation of the two-point density correlation along x- and z-direction for same $\phi = 0.44$ for different $p$ values. The minimum is slightly more pronounced for lower $p$ indicating smaller pore cluster in the z direction. This is highlighted in the inset.</figcaption></figure>

to lower values as porosity increases, as evident in Fig. 6(a). This shift is expected as the pore clusters have an elongated appearance at lower values of porosity (Figs 2a and b) and a greater separation between the pores. As porosity increases, the elongated pores start disappearing as greater correlation is introduced between the pores. The comparison of the density correlation along x- and z-directions in Fig. 6(a) clearly shows the anisotropy along the growth direction. It is expected that at $\mathbf{r} = 0$, the two-point density correlation must be equal to the porosity value corresponding to that $p$, as we have assigned a value of 1 to a pore site and a value 0 to a grain site. As $\mathbf{r} \to \infty$, here $\mathbf{r} = 200$ from Fig. 6(a), the two-point correlation function is expected to approach the value of $\phi^2$. This is also borne out from our calculations as is evident from Fig. 6(a).

Variation of the two-point correlation function with $\mathbf{r}$ calculated on the 3-D real rock sample in all the orthogonal directions is shown in Fig. 7. One unit of $\mathbf{r}$ is 5.06 µm in the real limestone sample. It seems that the real limestone sample is slightly anisotropic too. The bedding plane, referred to in the figure as the x-y plane, is isotropic in both directions but there is some anisotropy in the direction of deposition of grains. The direction of deposition of grains is the direction along which consecutive planes were cut and X-ray tomography taken. Since the two-point correlation for the bedding plane shows a lower minimum in the range of $\mathbf{r} = 0$ to 50, it indicates that the grain clusters in this plane have a larger size than along the growth direction. The density correlation at $\mathbf{r} = 0$ is 0.074 which matches the porosity of the limestone sample.

## 4 FRACTAL DIMENSION OF PORE VOLUME AND INTERFACE

Once the rock structure is generated, every cluster of pores is identified and checked for connectivity along the z-direction using the ‘burning algorithm’ of Hoshen & Kopelman (1976). The study of any transport property involves injection of fluid initially along the

growth direction as its perpendicular plane is exposed to the surface. Any transport can occur only through the connected pore cluster and so we have checked for connectivity only along this direction. But of course transport along other directions cannot be ruled out. A pore cluster is ‘connected’ if it has at least one pore on one surface and another pore of its cluster on the opposite surface of the sample. Only the six ‘nearest neighbours’ that meet the pore site along planes are checked for connectivity. In an earlier work of the authors on BBDM (Dutta & Tarafdar 2003), it was observed that only the ‘connected cluster’ was a volume fractal with the mass of the void space and the interface showing the same fractal dimension. When the calculation was done on the entire sample which included the disconnected pores also, the distribution of pores was found to be homogeneous with the same dimension as the Euclidean value. To investigate any impact of the compaction process, we checked for any signature of fractal nature in RBBDM for ‘total pore space’ and ‘connected pore space’ separately. The ‘pore centre’ was determined according to

$$ \mathbf{r}_p = \sum_{i=1}^{N} \frac{\mathbf{r}_i}{N}, \qquad (3) $$

where $\mathbf{r}_i$ is the position of the $i$th pore site and $N$ is the total number of pore sites. Here a vacant pixel site is counted as one pore site. The distribution of the pores when plotted against distance follows a power law

$$ N(r) \sim r^{d_f}, \qquad (4) $$

where $N(r)$ is the number of vacant cubes contained within a sphere of radius $r$ and centred about the ‘pore centre’ as obtained from eq.(3). $N(r)$ was calculated for concentric spheres of increasing radius $r$ measured with respect to the ‘pore centre’. Unlike BBDM, here the pore phase is a volume fractal when checked both for the entire structure inclusive of isolated pores and also for the ‘connected cluster’ phase. A log-log plot of the variation of $N_{\text{mass}}(r)$ with $r$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Two-point density correlation along three orthogonal directions of oolitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). A unit of *r* corresponds to 5.06 µm. The two almost coinciding graph are computed for direction parallel to bedding plane. The behaviour of red plot (perpendicular to bedding) near the origin is clearly different, suggesting a slight anisotropy of the medium.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Log–log plot of variation of pore mass and interface with distance from ‘void-centre’, shown for high porosity, $\phi_{max} = 0.457$ and low porosity $\phi = 0.231$. The inset shows the fractal dimension of pore mass and the interface for different values of *p*. It is clear that for all *p* values, the simulated rock is a volume fractal with the interface dimension equal to the pore mass dimension.</figcaption></figure>

for a high porosity value $\phi_{max} = 0.456$ and a low porosity value $\phi = 0.231$, for the connected pore cluster is depicted in Fig. 8. Real sedimentary rocks like sandstone and dolomite have an effective porosity between 0.21 and 0.27 (Pittmann 1984; Katz & Thompson 1985). This is higher than the porosity of the real limestone we have studied, but low compared to the $\phi_{max}$ as obtained from RBBDM. So we chose $\phi = 0.231$ as a case study of low porosity as shown in Fig. 8. The variation of the pore–solid interface $N_{int}(r)$ was also studied with increasing $r$. The slope of the log–log plot of $N_{int}(r)$ versus $r$ yielded the fractal dimension of the interface at any concentration $p$. In this case $N_{int}(r)$ is the number of cubes with the presence of any part of a interface contained within a sphere of radius $r$ measured with respect to the ‘pore centre’ of the spanning cluster. The log–log variation of $N_{int}(r)$ with $r$ for the same two values of porosities is

also shown in Fig. 8. It is evident that the graphs follow the pore mass curve almost exactly, the slope having the value between 2.8 and 2.95. This is greater than the Euclidian dimension of the interface, 2. This implies that the interface is extremely tortuous, almost volume filling, and has the same dimension as the mass as shown as an inset in Fig. 8. This is the characteristic of a volume fractal and is reminiscent of ‘peano curves’ (Vicsek 1992). The RBBDM is a volume fractal in the ‘total pore space’ and the ‘connected pore space’, unlike the BBDM. This difference with the BBDM is attributed to the compaction of the grains allowed during the growth process.

A closer examination of the graphs in Fig. 8 shows that each of these graphs show two regions with different slopes. This difference in slope is highlighted in Fig. 9 for the porosity 0.231 where the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Log–log plot of variation of pore mass with distance from ‘void-centre’, shown for porosity $\phi = 0.231$. The transition $r_{tr}$ is shown as the intersection of the two regression slopes shown by solid and dotted lines. The data points are shown as circles. The inset shows the two different fractal dimensions occurring in the two different length scales, plotted versus $p$. The difference decreases with higher $p$ values indicating that the sample becomes more homogeneous at lower porosities.</figcaption></figure>

pore mass is plotted versus $r$. The transition point $r_{tr}$, where the change in slope occurs, is determined as the intersection point of the two regression plots shown by solid and dotted lines in the graph. This change in slope, though small for every $p$ value studied, was distinct. In the inset to Fig. 9, we have plotted the two sets of fractal dimension of the pore mass obtained for each $p$ value. It is evident that the difference between the two sets of mass fractal dimension decreases with increasing $p$. The fractal dimension at distances greater than $r_{tr}$, all have values ~2.98 which is close to 3.0, the euclidean dimension. We conclude that the transition $r_{tr}$ indicates the distance at which the structure for each $p$ changes over from its fractal nature to euclidean nature. From Fig. 2(a) we see that though the pore clusters show more homogeneity, they retain an elongated structure along the direction of assembly of grains whereas from Figs 3(a) and (b), we see that the pore clusters are almost channel-like with practically no connectivity in the direction perpendicular to direction of assembly. With channel-like pores, the pore space is expected to have a euclidean nature. In terms of porosity, we conclude that the simulated rock sample shows greater inhomogeneity in pore shape at higher porosity (lower $p$ value) and tend to become more homogeneous at lower porosity (higher $p$). It seems that the sedimentary rock simulated by RBBDM may have multifractal characteristic.

The variation of pore mass and interface with distance on the reconstructed real rock structure is depicted in Fig. 10. Though the variation of pore mass and interface with distance from the ‘void centre’, closely follow one another on a log–log plot, a clear linear region is not discernible. Real rocks may not be a monofractal, and this possibility has to be investigated further. However the fact that the pore mass and interface follow each other so closely, may be indicative that the pore space in the real rock studied, is a volume fractal!

The measurements on porous rocks usually report a fractal dimension for the pore–solid interface rather than the pore volume. Adsorption experiments essentially look only at the interface. The measured area of the interface depends on the length scale used as probe, over several orders of magnitude of probe size. Scattering

experiments which can identify volume fractals as well as surface fractals, usually indicate surface fractal results for rocks. The neutrons or X-ray photons are scattered by the solid grains, so the probe in this case sees the solid grains. This information should be interpreted as follows: the pore–solid interface is fractal, the pore volume may or may not be fractal while the solid rock phase is non-fractal (Sen et al. 1981; Avnir et al. 1984; Bale and Schmidt 1984; Schaefer and Keefer 1984; Wong et al. 1984; Katz and Thompson 1985; Sen et al. 2002).

## 5 ANOMALOUS DIFFUSION IN THE PORE SPACE

To study diffusion, we follow the ‘blind ant’ algorithm (Havlin & Ben-Avraham 1987) on the connected cluster, also referred to as the ‘infinite cluster’ in random walk problems. It makes sense to restrict ourselves to infinite clusters, as the study of diffusion on finite clusters may result in walkers spending infinite time trapped in such clusters, giving rise to misleading results in the study of diffusion. The ‘blind ant’ starts its random walk from a vacant site chosen randomly from the mid-section of the generated rock sample to ensure that the ant was allowed sufficient time before it could reach any of the boundaries and terminate its walk. Periodic boundary conditions are not introduced as this would destroy any self-similarity of the fractal space. Every configuration was allowed 500 independent walks before the next new configuration was chosen. The ant chooses any one of its nearest neighbours with equal probability and attempts to jump to that site. If the chosen site is not available (in case it is a solid grain site), it waits in its original site and again attempts to make a jump in the next time step. So the blind ant may or may not move in every time step. The variation of the mean square distance $\langle r^2 \rangle$ traversed in a given time $t$ gives the diffusion equation. The displacement of the ant is measured from its starting position and every attempt to jump is an increase in time. The mean square displacement follows a relation

$$ \langle r^2 \rangle = \alpha t^\nu, \qquad (5) $$
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Log-log plot of variation of pore mass and interface with distance from 'void-centre' of oolitic limestone(pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). A unit of *r* corresponds to 5.06 µm. The two plots follow each other, suggesting that the pore space of the oolitic limestone is possibly a volume fractal.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Variation of mean square displacement along *x*- and *z*-directions, with time on connected pore cluster of simulated rock structure at low and high porosity. This displacement is greater in the *z*-direction as the pore clusters are more elongated in that direction. The anisotropy decreases as the porosity increases.</figcaption></figure>

where **r** is the displacement in time *t*. For normal diffusion, *v* = 1 and the proportionality constant *α* would be the diffusion coefficient. For *v* < 1, the diffusion is anomalous and this is expected in fractal structures (Havlin & Ben-Avraham 1987). Mean square displacement along the *x*- and *z*-axes of the simulated structure for a low and a high porosity value, is plotted in Fig. 11. The mean square displacement along the *y*-direction was similar to the mean square displacement along the *x*-direction. This once again shows that though the simulated structure shows isotropy in the transverse plane, it is anisotropic in the direction of growth. This anisotropy decreases as the porosity increases. A log-log plot of ⟨r²⟩ versus *t* for different *p* values is shown in Fig. 12. The slope of each plot gives the corresponding *v* value. The *v* values show variation with

*p*, though they all have values <1 indicating anomalous diffusion. This is expected in a fractal pore space.

The variation of ⟨r²⟩ versus *t* for a high porosity *φ* = 0.45 is depicted in Fig. 13 to highlight that there is a change of slope at time *t*<sub>r</sub>, which is termed transition time. It is difficult to relate the variation of *t*<sub>r</sub> to the microstructure at this stage. It is possible that the structure could have a multifractal nature which shall be investigated in the future. The *v* values before and after *t*<sub>r</sub> for any *p*, are plotted against *p* in the inset of Fig. 13. The *v* values before and after *t*<sub>r</sub> show two clear clusters which indicate the presence of multifractal nature. The difference between the *v* values is more pronounced at lower *p* values where the porosity is high. As *p* increases the difference decreases and approaches 0.85 at *p* = 0.9, indicating that
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Log-log plot of mean square displacement with time on connected pore cluster of simulated rock structure for different *p*.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Log-log plot of mean square displacement with time on connected pore cluster of simulated rock structure to highlight the existence of $t_r$. The linear regions on either side of $t_r$ are shown by solid and dotted lines. $t_r$ is evaluated from the intersection of the two slopes. The data points are shown as circles. The inset shows the variation of the exponent $v$ with $p$. $v$ for all $p$ is always <1, indicating anomalous diffusion. There are two clear clusters of $v$ below and above $t_r$.</figcaption></figure>

the structure tends to be a monofractal here. Between $p = 0.9$ and 1.0, $v$ approaches 1, implying that diffusion tends to be normal. With reference to Fig. 3 we may conclude that since at very low porosities corresponding to $p > 0.9$, the pore clusters are straight channel-like, diffusion through them is expected to be normal.

Diffusion was studied on the 3-D real oolitic limestone(pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). This was generated from the array of 2-D sections of X-ray microtomographs. The variation of mean square displacement with time along x- and z-direction is plotted in Fig. 14. The mean square displacement along the y-direction behaviour is almost the same as along the x-direction. The real rock sample also shows some amount of anisotropy, which may be attributed to the difference in the micro geometry of the pore clusters. This is corroborated by

our study of the two-point correlation, refer to Fig. 7, and discussed in the section on it. The variation of $\langle r^2 \rangle$ versus $t$ on the real rock sample is shown in Fig. 14. Two different slopes are evident with $v = 0.926$ in the first region and $v = 0.665$ in the second region. This difference remained even though several walks were carried out on this structure with the walk starting from different positions within the sample. Normal diffusion shows $v = 1.0$ which is close to $v = 0.926$ observed within a mean square displacement $\sim 10$ units. This indicates that diffusion is closer to normal here. An examination of Figs 4(a) and (b) clearly shows large white regions, which are the large voids between the black grains whose outlines are visible. The white spaces are multiply connected clusters of vacant pixels, which have a considerable volume, akin to the 'pore bodies' in earlier models (Sen *et al.* 1981). In the present analysis of diffusion,
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Log-log plot of mean square displacement with time along x- and z-orthogonal directions, on connected pore cluster of oolitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). A unit of *r* corresponds to 5.06 µm. The anisotropic nature of the pore space is clear from the similarity of the plots. Log-log plot of mean square displacement versus time on the oolitic limestone shows two linear regions with *v* values 0.926 and 0.665.</figcaption></figure>

a random walker within such a pore body will experience normal diffusion with *v* close to 1. Given a time which exceeds the time required to explore the pore body, the walker has to move through narrow channels, barely visible in Figs 4(a) and (b), to reach another large pore. This time is a ‘crossover time’, after which the walker ‘sees’ a complex tortuous and non-euclidean or fractal space. This explains the two different exponents *v*, one corresponding to normal diffusion in euclidean space and the other to anomalous diffusion in a fractal space.

An estimate of the average pore cluster size may be done from diffusion study on the real rock sample. The mean square displacement is related to diffusion time by eq. (5). Transition time *t*, corresponds to ⟨r²⟩ ~ 10 from Fig. 14. Therefore, since

$$⟨r^2⟩ = 10$$

$$⇒ ⟨r⟩ = 3.2 \text{ pixels (1 unit} \equiv \text{ 1 pixel)}$$

$$⇒ ⟨r⟩ \sim 16 \mu\text{m}$$

Thus the average pore body ⟨r<sub>body</sub>⟩ is ~ 16 µm. This matches the average value of 4 pixels, that is, ~20 µm obtained from processed X-ray microtomographs, thus justifying our analysis.

Beyond this distance, the average pore size distribution is heterogeneous and fractal in nature.

## 6 DISCUSSION

We have investigated the pore space of a sedimentary rock using the RBBDM for different values of porosity and compared our results with a real system. The RBBDM is an improvement on the BBDM as it includes the compaction of the grains during the deposition process. The compaction is mimicked through the toppling of unstable overhangs that develop during growth. Compared to our previous model BBDM, it is evident that compaction has resulted in a non-monotonic change of porosity with *p*, the fraction of smaller grains. In this case, the maximum porosity was obtained at *p* = 0.5. In an earlier study (Dutta & Tarafdar 2005), the authors had generated a rock structure with RBBDM using an aspect ratio of

1: 3 for the larger grain. No discernible difference in results was noted. However with a different combination of grain size, it may be possible to have different variation of porosity values as the correlation that is introduced between adjacent columns of grains will be affected differently through the toppling process.

Another important effect of compaction in our model is that the both the ‘connected’ and the ‘total’ pore space of the sedimentary rock is a volume fractal. In the BBDM, only the ‘connected’ pore space was found to be a volume fractal, whereas the total pore space was Euclidean in nature. The process of toppling of the grains ensures that the interface always remains highly convoluted, almost space filling, that is it has the same dimension as the embedding space. A convenient way of distinguishing a tortuous space filling fractal from a normal Euclidean structure, is to look at the diffusion of a random walker through the structure. Diffusion will be anomalous in the latter case. This is exactly what we find here.

The exponent *v* has a value <1 for all *p* values, indicating that diffusion is subdiffusive. This was also observed in our study of diffusion in the pore space of the reconstructed real rock obtained from compiling 1000 consecutive 2-D microtomographs. So RBBDM is a good model for generating sedimentary rocks of different porosities as it produces similar features of pore space. The study of two-point density correlation of the pores on the RBBDM shows the same nature as obtained in the real rock. However the growth rule of the RBBDM results in an anisotropy in the pore space, the pore clusters are more elongated in the direction of growth than in its transverse plane. This anisotropy is manifested in the density correlation study and also in the component wise mean square displacement in diffusion. The anisotropy decreases as the fraction of smaller grains, that is, *p* value decreases. The real rock sample was also found to show anisotropy in correlation and diffusion studies. The bedding plane was found to be isotropic in the two perpendicular directions, but exhibited anisotropy in the direction of deposition.

Investigation of the pore mass and pore-rock interface indicated that sedimentary rocks may have multiscaling properties. Studies on both the simulated and real rock structures suggested this possibility.
Our diffusion studies also strengthened this possibility. Whether the sedimentary rock has a multifractal nature, has to be investigated thoroughly before a clear conclusion can be drawn.

An estimate of average pore size was calculated from diffusion study through the real rock studied. The result matched the data obtained from X-ray tomography.

Our studies show that the RBBDM is a good model for simulating sedimentary rocks as its characteristics are similar to real rock samples. The characterization of the pore space of the RBBDM will help us to correlate the transport properties of any fluid through it, to its microgeometry. Our future plans are oriented in this direction.

## ACKNOWLEDGMENTS

This work is supported by Indo-French Centre For the Promotion of Advanced Research (IFCPAR project no:4409-1). AG is grateful to IFCPAR for providing a research fellowship.

## REFERENCES

* Avnir, D., Farin, D. & Pfeifer, P., 1984. Molecular fractal surfaces, *Nature*, **308**, 261–263.

* Bale, H.D. & Schmidt, P.W., 1984. Small angle X-ray scattering of submicroscopic porosity with fractal properties, *Phys. Rev. Lett.*, **53**, 596.

* Berryman, J.G. & Blair, S.C., 1986. Use of digital image analysis to estimate fluid permeability of porous materials: application of two-point correlation functions, *J. appl. Phys.*, **60**, 1930.

* Dutta, T. & Tarafdar, S., 2003. Fractal pore structure of sedimentary rocks: simulation by ballistic deposition, *J. geophys. Res.*, **108**(B2), 2062.

* Dutta, T. & Tarafdar, S., 2005. Simulation of sedimentary rocks by ballistic deposition, in *Proceedings of 9th World Multiconf. on Systemics, Cybernetics and Informatics*, Orlando, FL, USA, 10–13, p. 151.

* Havlin, S. & Ben-Avraham, D., 1987. Diffusion in disordered media, *Adv. Phys.*, **36**(6), 695.

* Hornung, U., 1997. *Homogenization and Porous Media*, Springer-Verlag Inc, New York.

* Hoshen, J. & Kopelman, R., 1976. Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm, *Phys. Rev. B*, **14**, 3438.

* Karmakar, R., Dutta, T., Lebovka, N. & Tarafdar, S., 2005. Effect of surface roughness on the bulk properties of simulated porous media, *Physica A*, **348**, 236–244.

* Katz, A.J. & Thompson, A.H., 1985. Fractal sandstones pores: implications for conductivity and pore formations, *Phys. Rev. Lett.*, **54**, 1325.

* Krohn, C.E., 1988a. Fractal measurements of sandstones, shales, and carbonates, *J. geophys. Res.*, **93**, 3297.

* Krohn, C.E., 1988b. Sandstone fractal and Euclidean pore volume distributions, *J. geophys. Res.*, **93**, 3286.

* Manna, S.S., Dutta, T., Karmakar, R. & Tarafdar, S., 2002. A percolation model for diagenesis, *Int. J. Mod. Phys. C*, **13**, 319.

* Manwart, C., Torquato, S. & Hilfer, R., 2000. Stochastic reconstruction of sandstones, *Phys. Rev. E.*, **62**, 893.

* Pettijohn, F.J., 1984. *Sedimentary Rocks*, Harper & Row Publishers Inc., USA.

* Pittmann, E.D., 1984. *Physics and Chemistry of Porous Media*, AIP, New York, p. 1.

* Sadhukhan, S., Dutta, T. & Tarafdar, S., 2007a. Simulation of diagenesis and permeability variation in two-dimensional rock structure, *Geophys. J. Int.*, **169**, 1366–1375.

* Sadhukhan, S., Dutta, T. & Tarafdar, S., 2007b. Pore structure and conductivity modelled by bidisperse ballistic deposition with relaxation, *Modeling Simul. Mater. Sci. Eng.*, **15**, 773–786.

* Sadhukhan, S., Mal, D., Dutta, T. & Tarafdar, S., 2008. Permeability variation with fracture dissolution: role of diffusion vs. drift., *Physica A*, **387**, 4541–4546.

* Schaefer, D.W. & Keefer, K.D., 1984. Fractal geometry of silica condensation polymers., *Phys. Rev. Lett.*, **53**, 1383.

* Sen, D., Mazumder, S. & Tarafdar, S., 2002. Pore morphology and pore surface roughening in rocks: a small angle neutron scattering investigation, *J. Mater. Sci.*, **37**, 941.

* Sen, P.N., Scala, C. & Cohen, M.H., 1981. A self-similar model for sedimentary rocks: application to the dielectric constant of fused glass beads, *Geophys.*, **46**, 781.

* Stauffer, D. & Aharony, A., 1994. *Introduction to Percolation Theory*, 2nd edn, Taylor and Francis London, UK.

* Tarafdar, S. & Roy, S., 1998. A growth model for porous sedimentary rocks, *Physica B*, **254**, 28.

* Vicsek Tamas, 1992. *Fractal Growth Phenomena*, 2nd edn, World Scientific Publishing Co. Pvt. Ltd Singapore.

* Whitaker, S., 1999. *The Method of Volume Averaging*, Kluwer Academic Publishers the Netherlands.

* Wong, P.Z., Koplik, J. & Tomanic, J.P., 1984. Conductivity and permeability of rocks, *Phys. Rev. B.*, **30**, 660.

* Wong, P.Z., 1987. Fractal surfaces in porous media, in *Proceedings of the 2nd International Symposium on Physics and Chemistry of Porous Media, AIP Conf. Proc.*, **154**, 304.