# Incorporating metamorphism in geodynamic models: the mass conservation problem

Gyoergy Hetenyi, Vincent Godard, Rodolphe Cattin, James A. D. Connolly

► To cite this version:

Gyoergy Hetenyi, Vincent Godard, Rodolphe Cattin, James A. D. Connolly. Incorporating metamorphism in geodynamic models: the mass conservation problem. Geophysical Journal International, 2011, 186 (3), pp.6-10. 10.1111/j.1365-246X.2011.05052.x. hal-00617480

HAL Id: hal-00617480

https://hal.science/hal-00617480v1

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 CC BY 4.0 - Attribution - International License
# Multifractal analysis of the pore space of real and simulated sedimentary rocks

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 2014 October 22. Received 2014 October 21; in original form 2014 May 28

## SUMMARY

It is well known that sedimentary rocks having same porosity can have very different pore size distribution. The pore distribution determines many characteristics of the rock, among which its transport properties are often the most useful. Multifractal analysis is a powerful tool that is increasingly used to characterize the pore space. In this study, we performed multifractal analysis of pore distribution on simulated sedimentary rocks using the relaxed bidisperse ballistic deposition model (RBBDM). The RBBDM can generate a 3-D structure of sedimentary rocks of variable porosity by tuning the fraction *p* of particles of two different sizes. We also performed multifractal analysis on two samples of real sedimentary rock to compare with the simulation studies. One sample, an oolitic limestone is of high porosity (40 per cent) while the other is a reefal carbonate of low porosity, around 7 per cent. 2-D sections of X-ray microtomographs of the real rocks were stacked sequentially to reconstruct the real rock specimens. Both samples show multifractal character. The results from analysis of real rock agree quite well with the simulated structure of low porosity. The simulated rock of high porosity showed a weak multifractal nature though the real rock sample of similar porosity was found to be strongly multifractal. We propose a ‘structure parameter’ ξ which is a function of porosity and the generalized dimensions, and controls the transport properties of the rock.

**Key words:** Spatial analysis; Fractals and multifractals; Microstructures.

## 1 INTRODUCTION

Sedimentary rocks are often the storehouses of natural oil and gases whose extraction depend on the permeability of these fluids through them. The transport properties of sedimentary rocks depend not only on the porosity of the rocks but more importantly on the pore size distribution (PSD) and their connectivity. The pore space can be a continuum of pores with extremely varying pore sizes ranging over a scale of $10^6$, besides being extremely complex and heterogeneous and often self-similar.

Fractal and multifractal analysis are increasingly used to study complex heterogeneous systems which show self-similarity on several length scales. They have the ability to provide an accurate representation of the heterogeneous pore geometry and address the relationship between porosity and a range of physical processes happening in a porous medium like transport of water in soils, extraction of oil and natural gases and CO₂ sequestration in sedimentary rock. The role of structure of the pores and their size distribution may result in fractal structures like viscous fingering during biphasic flow occurring in porous media (Toussaint et al. 2005; Lovoll et al. 2011; Nsir et al. 2012). Multifractal analysis

on porous media has been done using fractal models (Hansen & Skjeltorp 1988; Hansen et al. 1988; Rieu & Sposito 1991), image analysis of 2-D sections of soil blocks (Tarquis et al. 2003; Dathe et al. 2006; Grau et al. 2006), analysis of 3-D pore systems reconstructed by computer tomography (Tarquis et al. 2007), mercury intrusion porosimetry (Vidal Vázquez et al. 2008) and nitrogen absorption isotherms (Paz-Ferreiro et al. 2009).

In this work, the authors use the relaxed bidisperse ballistic deposition (RBBDM) to simulate a 3-D porous rock structure of varying porosity. In an earlier effort (Giri et al. 2012) to probe the geometry of the microstructure of the pore clusters produced by the RBBDM at different porosities, the authors had noticed that the simulated structure had a fractal nature that changed across a characteristic length $r_{fr}$. The power-law exponent of the pore mass had different values across $r_{fr}$, hinting that the pore space might have a multifractal nature. This idea was further strengthened by diffusion studies in connected pore clusters. For the entire range of porosities studied, diffusion was found to be anomalous with different values of diffusion exponent over different length scales. Real rock samples were studied for comparison of simulation results, and similar signature of multifractal nature was found there.

<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>
In this study, our objectives are: (i) to investigate whether the simulated porous structure generated with the RBBDM is indeed a multifractal in its pore geometry; (ii) to bring out differences in the microgeometry of the simulated structure at different porosities through a study of their multifractal spectral dimensions; (iii) to relate the multifractal parameters at different length scales to macroscopic properties like transport properties of the rock; (iv) to compare our results with similar studies done on real limestone and carbonate rock samples and (v) to establish a link between the multifractal characteristics and transport properties in rocks.

# 2 MATERIALS AND METHODS

## 2.1 Sedimentary rock samples

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'. The details of RBBDM have been discussed in earlier works (Sadhukhan et al. 2007a,b, 2008) in order to study various transport properties like permeability and conductivity through sedimentary rocks. A brief outline of the model will be given here. The basic algorithm is to deposit particles of two different sizes ballistically. In three dimensions (2 + 1 model), we drop square $1 \times 1 \times 1$ and elongated $2 \times 1 \times 1$ grains on a square substrate. The longer grains deposit with their long axis parallel to the bedding plane and along either of the two tranverse directions with equal probability. Natural sand grains are angular and elongated (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 cubic grains are chosen with a probability $p$ and elongated grains with probability $(1-p)$. The porosity $\phi$, defined as the vacant fraction of the total volume, depends on the value of $p$. For $p=1$, a dense structure with no pores is produced. When $p$ is gradually decreased the presence of a single large grain sitting atop a long pore cluster, introduces correlation between adjacent columns (Karmakar et al. 2005). So a substrate of sufficient height needs to be generated before the porosity value can stabilize (Manna et al. 2002; Dutta & Tarafdar 2003). 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 if there is no supporting particle immediately below the protrusion of the second overhang. A two-step overhang created in this way, is not stable and the second large particle topples over if possible, according to the rule scheme as shown in Fig. 1 leading to compaction.

<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.** (a) $x-z$ section of simulated structure for $\phi_{max} = 0.45$, that is high porosity. (b) $x-y$ section at same porosity. Structure looks more isotropic. The white indicate pore clusters.</figcaption></figure>

In their earlier works (Manna et al. 2002; Dutta & Tarafdar 2003), the authors have discussed that the sample attains a constant porosity only when sufficient number of grains (depending on sample size) has been deposited. Here a $Lx \times Ly \times Lz$ size sample was generated, from which a $Lx \times Ly \times Lz$ sample was selected after the porosity had stabilized to within 0.001 per cent. The selected sample was chosen from below the deepest trough at the surface to eliminate surface effects. All results on simulation are reported for $512 \times 512 \times 512$.

A vertical section, $(x-z)$plane, of the generated sample at porosity $\phi = 0.42$ is shown in Fig. 2(a), while a horizontal section $(x-y)$ plane) of the sample at the same porosity value is shown in Fig. 2(b). The anisotropy in the pore geometry is clearly visible. The pore clusters have an elongated and interconnected appearance along the $z$-direction while the distribution of pores along the horizontal plane is quite homogeneous. As the fraction of larger grains is decreased, the porosity of the sample decreases and the pore distribution becomes more anisotropic nature. The vertical $(x-z)$ and $(y-z)$ planes, and the horizontal bedding section, $(x-y)$ plane, of the sample at a very low porosity $\phi = 0.077$, are shown, respectively in Figs 3(a) and (b) and (c). The elongated isolated pore clusters are prominent in the direction of assembly of the grains, whereas the pores remain homogeneously distributed in the $(x-y)$ plane.

To compare our simulation results with real rock samples, X-ray tomography micrographs of 2-D sections of two real sedimentary rock samples obtained from an oolitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France), and a reefal carbonate from the Majorca Island, Spain, have been used. The oolitic limestone is composed of recrystallized ooliths with a mean diameter of less than a few hundred µm. Each pixel of both the micrographs corresponds to 5.06 µm. For every real rock sample studied, each micrograph section was converted to a binary file form such that 0 corresponded to a pore site and 1 corresponded to a rock site. An array of 1000 consecutive binary files were put together precisely to reconstruct the form of the real 3-D rock structures. To obtain a picture of the real rock sample, each binary file was converted to a grey-scaled picture, as shown in Fig. 4, using ImageJ. In each of the two samples, the real structure chosen was $1000 \times 1000 \times 1000$ in size, and all studies on real rock was carried out on this structure. The 2-D sections of oolitic limestone cut in the direction of assembly (growth) from the reconstructed 3-D structure are shown in Fig. 4(a) [(x-z) plane] and Fig. 4(b) [(y-z) plane]. Comparison between Figs 4(a) and (b) shows the pore distribution to be isotropic. A 2-D section of the bedding plane of the same rock structure is shown in Fig. 4(c).
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) and (b) Section of simulated structure along x-z planes and y-z planes, respectively for $\phi = 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. (c) A section of x-y plane at same porosity. Anisotropy in pore cluster structure is quite pronounced. The white indicate pore clusters.</figcaption></figure>

A slight anisotropy between the pore distribution in the bedding plane and the direction of growth is discernible from a comparison between Figs 4(a)–(c). Similar 2-D sections of the reefal carbonate rock sample have been cut along the direction of assembly of the reconstructed three-dimension sample and shown in Figs 4(d) and (e) while a section of the same sample cut along the bedding plane is displayed in Fig. 4(f). Unlike the previous sample, any anisotropy that may be present in the PSD along the two directions, is not easily discernible.

## 2.2 Multifractal concepts

Highly inhomogeneous systems which do not obey self-similar scaling law with a single exponent, may actually consist of several intertwined fractal sets with a spectrum of fractal dimensions. These systems are said to be multifractal. Such systems have a complex distribution which arises from peculiarities of their generation. These are not a simple collection of fractal systems, rather one may say that these constitute a distribution of several fractal subsets on a geometrical support. Each of these subsets is characterized by its singularity strength and fractal dimension.

Multifractal analysis involves the estimation of three functions: mass exponent $\tau_q$, singularity strength (or local scaling index) $\alpha_q$, and multifractal or singularity spectrum $f(\alpha)$. We follow the procedure for multifractal analysis as elaborated in standard texts on fractals (Feder 1988; Vicsek 1992).

The chosen system of size $L^3$ is divided into a set of different boxes (cubes) of equal size having length $l$. To make the analysis system size independent, we define a characteristic length $\epsilon = l/L$. For a system of size $512 \times 512 \times 512$, a common choice is to consider dyadic scaling down, that is successive partitions of the system in $k$ stages ($k = 1, 2, 3, ...$), that generate a number of cells $N(l) = 2^{3k}$ of characteristic length $l = L \times 2^{-k}$. For images of size $512 \times 512 \times 512$ pixels, choosing $\epsilon$ in this manner avoids

artefacts which occur when boxes do not entirely cover the image at the borders. In terms of $\epsilon$, $N(\epsilon)$ are the number of boxes of characteristic length $\epsilon$ required to cover the entire system. Each such box is labelled $i$.

The probability mass function $\mu_i(\epsilon)$ describing the portion of the ‘mass’ contained in the $i$th box of size $\epsilon$ is given by

$$ \mu_i(\epsilon) = \frac{m_i}{m_i}, \qquad (1) $$

where $m_i$ is the number of pore sites in the $i$th box and $m_i$ is the total number of pore sites in the entire system. A pore site refers to a pixel that is vacant.

The partition function $\chi(q, \epsilon)$ for different moments $q$ is estimated from $\mu_i^q(\epsilon)$ values as

$$ \chi(q, \epsilon) = \sum_{i=1}^{N(i)} \mu_i^q(\epsilon). \qquad (2) $$

The box size $\epsilon$ may be considered as a filter so that by changing $\epsilon$ one may explore the sample at different scales. The partition function $\chi(q, \epsilon)$ contains information at different scales and moments. The sum in the numerator is dominated by the highest value of $\mu_i$ for $q > 0$ and the lowest value of $\mu_i$ for $q < 0$.

The measure of the $q$th moment of the mass distribution of the system is defined as

$$ M(q, \epsilon) = \sum_{i=1}^{N} \mu_i(\epsilon)^q \epsilon^d. \qquad (3) $$

If $\sum_{i=1}^N \mu_i^q(\epsilon)$ in the limit $\epsilon \to 0$ crosses over from 0 to $\infty$ as $d$ changes from a value less than $\tau(q)$ to a value greater than $\tau(q)$, then the measure has a mass exponent

$$ d = \tau(q). \qquad (4) $$

The dimension of the $q$th exponent of the mass distribution is obtained by counting the number of boxes $N(q, \epsilon)$ of size $\epsilon$ needed to cover the 3-D structure whose multifractal character is being determined. Using eq. (2), partition function $\chi(q, \epsilon)$ for different moments $q$ is estimated. The partition function $\chi(q, \epsilon)$ scales on varying $\epsilon$ as

$$ \chi(q, \epsilon) = \epsilon^{-\tau(q)}. \qquad (5) $$

So the mass exponents $\tau(q)$ of the measure for different $q$ are determined. Here $q$ can vary from $+\infty$ to $-\infty$.

The $q$th moment of the probability mass function as defined in eq. (1) may be represented as $\mu_i(q, \epsilon)$, and it also scales with $\epsilon$ (Halsey et al. 1986; Chhabra et al. 1989) as

$$ \mu_i(q, \epsilon) = \epsilon^{\alpha_i(q)}, \qquad (6) $$

where $\alpha_i(q)$ is the Hölder exponent or ‘singularity exponent’ or ‘crowding index’ of $\mu_i(q, \epsilon)$ particular to each $i$th box. A higher value of Hölder exponent implies a smaller concentration, and vice versa. The weighted singularity strength with respect to $\mu_i(q, \epsilon)$ can be obtained as an average of the $\alpha_i$’s weighted by $\mu_i(q, \epsilon)$ evaluated in the limit

$$ \alpha(q) = \lim_{\epsilon \to 0} \frac{\sum_{i=1}^{N(\epsilon)} \mu_i(q, \epsilon) \log \mu_i(q, \epsilon)}{\log \epsilon} \qquad (7) $$

whenever this limit exists.

The Hölder exponents $\alpha$ becomes a non-increasing function of $q$. Large (small) values of the parameter $q$ correspond to high (low)
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Each section is a square of side 2.58 × 10⁻³ m. The black and the white indicate rock and pore phases. (a) and (b) Section of the real limestone cut in mutually perpendicular planes. These are perpendicular to bedding planes. (c) A section of the bedding plane. These are sections of olitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France). (d) and (e) Section of reefal carbonate obtained from Majorca Islands, Spain, perpendicular to bedding planes. (f) A section of the carbonate rock along the bedding plane. This rock structure looks more isotropic.</figcaption></figure>

degrees of concentration of the measure. This is a natural result: α(q) is obtained as an average with respect to the probability measure μ(q) that magnifies the denser (more rarefied) regions for large (small) values of q. In particular, α(0) is the average of the Hölder exponents weighted by the uniform distribution over the support of μ; so that the sets where the measure μ is most concentrated are weighted like those where it is most rarefied; however, the average is weighted by the probability measure μ itself in the case of α(1). The analysis of the value of α(q) for different values of the parameter q will serve to evaluate the singular nature of the experimental measure μ: the larger variation of α(q) with respect to q; the greater singularity degree of the measure μ. Singularity exponents of multifractal distributions show a great variability within an interval (α_min, α_max) when ϵ tends to zero. For a monofractal, this interval reduces to a point.

Again, the number N_ϵ(α) of boxes of size ϵ that have a Hölder exponent between α and α + δα obeys a power law as

$$N_{\epsilon}(\alpha) \propto \epsilon^{-f(\alpha)}, \qquad (8)$$

where f(α) is a scaling exponent of the boxes with a common α, called the singularity exponent. A plot of f(α) versus α is called the singularity spectrum. f(α) is the fractal dimension of the set of points that have the same singularity exponent α. There can be several such interwoven fractal sets of points each with its particular value of f(α). Within each such set, the measure shows a particular scaling described by α.

Following (Halsey et al. 1986), the functions α and f(α) can be determined by Legendre transformation as

$$\alpha(q) = - \frac{d\tau(q)}{dq} \quad \text{and} \quad f(\alpha) = \alpha(q)q - \tau(q). \qquad (9)$$

Thus the singularity exponent defined by eq. (6) becomes a decreasing function of q. As q varies, points [α_q, f(α_q)] define a parabolic curve that attains a maximum value f(α_0) at the point α_0. f(α_0) gives the fractal dimension of the support as obtained by the box-counting method (Vicsek 1992).

Another equivalent description of the multifractal system is obtained from D_q versus q plot, where D_q, called the generalized dimension, corresponds to the scaling exponent for the qth moment of the measure. It is defined (Rényi 1960) by

$$D_{q \to 0} = \frac{1}{1-q} \frac{\log[\chi(q, \epsilon)]}{\log(\epsilon)}. \qquad (10)$$

For the particular case of q = 1, eq. (10) becomes indeterminate, and is estimated by l'Hôpital's rule

$$D_1 = \lim_{\epsilon \to 0} \frac{\sum_i \mu_i \ln \mu_i}{\ln \epsilon} \qquad (11)$$

D_q is related to the mass exponent τ(q) by

$$\tau(q) = (1 - q)D_q. \qquad (12)$$

In general, D_q is a non-increasing function that is constant in the case of standard smooth distributions. The higher the degree of heterogeneity of the measure, the greater is the variation in the values of D_q. For q ≫ 1, the dimension quantifies the scaling
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Pore size distribution for real and simulated rock structures. (a) Simulated rock: low porosity $\phi \simeq 0.07$. (b) Simulated rock: high porosity $\phi \simeq 0.42$. The inset shows the maximum contribution to PSD shown in main figure. (c) Limestone rock: low porosity $\phi \simeq 0.07$. (d) Carbonate rock: high porosity $\phi \simeq 0.4$.</figcaption></figure>

behaviour of denser cells, while for $q \ll 1$, the dimension measures the scaling properties of the rarer ones. The generalized dimensions $D_q$ for $q = 0, q = 1$ and $q = 2$ are known as the capacity, the information (Shannon entropy; Shannon 1949) and correlation dimensions (Grassberger & Procaccia 1983), respectively. Mathematically, the multifractals can be completely determined only by the entire multifractal spectrum. The values $D_q$ thus provide a valuable characterization of the singular behaviour of the measure and the respective interpretation within each context. However we use a few characteristic functions to describe the main properties of multifractals.

## 3 RESULTS AND DISCUSSION

The RBBDM has the potential for generating a porous three-dimension rock structure of porosity raging between 0 and 1. To test for any multifractal character that a simulated sample may have, we could present the analysis on rocks having any porosity value within this range. Though we have checked for possible multifractal signature of porous rocks having different porosity values by changing $p$, in this paper we present our detailed analysis on only two simulated rock samples, one of low porosity value, $\phi = 0.07$, and another of high porosity, $\phi = 0.42$. We deliberately choose these so as to compare the results with the two real rocks samples we had, which had similar porosity values.

### 3.1 Multifractal analysis of simulated rock structure

In the case of the simulated structure, a 512 x 512 x 512 cube was selected from below the deepest trough from the surface after the porosity had stabilized in an initial structure of 512 x 512 x 5000,

and after the porosity had stabilized for a particular choice of $p$. We calculated the size of all pore bodies present in the structure generated and plotted the number of separate pore bodies of a particular size versus the frequency of its occurrence to obtain the PSD. A comparison of the PSDs of the simulated structure for a low porosity value and a high porosity value, shown in Figs 5(a) and (b), respectively, show a larger size distribution of pore bodies in the rock with a lower porosity than that of high porosity. The narrow distribution of pore sizes for the high porosity sample is magnified in the inset of Fig. 5(b).

To carry out the multifractal analysis the system was covered with cubes of length size $\epsilon = 2^k$ with $k$ ranging between 1 and 9. The partition function $\chi(q, \epsilon)$, calculated according to eq. (2) for different values of box size $\epsilon$ and for different moment values $q$ was determined. A log-log plot of $\chi(q, \epsilon)$ versus $\epsilon$ when plotted for both low and high porosity samples, Figs 6(a) and (b), respectively, showed a deviation from linearity beyond a certain range of $q$. A power law scaling was observed in the range $q = -9$ to +9. All calculations have been done within this range of $q$ where the coefficient of determination ($R^2$) was greater than 0.98 for both the simulated samples as shown in Table 1. The exponent $\tau(q)$ for each such $q$ and for both samples was noted.

The scaling properties observed in the partition function can be characterized by determining if the scaling is simple as in monofractal, or multiple as in multifractal. The variation of $\tau(q)$ versus $q$ for a low porosity $\phi = 0.08$ corresponding to $p = 0.997$, and a high porosity $\phi = 0.42$ corresponding to $p = 0.8$ are shown in Figs 7(a) and (b), respectively. The data points from simulation studies are shown as open circles in the graphs. The $\tau(q)$ functions which are straight lines for monofractals, deviate from linear behaviour for multifractals. The slopes of the $\tau(q)$ for $q < 0$ are different from those for $q > 0$, the difference being greater for low porosity than
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Log-log plot of $\chi$ versus $\epsilon$ for different moments. (a) Simulated rock: low porosity $\phi \simeq 0.07$. (b) Simulated rock: high porosity $\phi \simeq 0.42$. (c) Limestone rock: low porosity $\phi \simeq 0.07$. (d) Carbonate rock: high porosity $\phi \simeq 0.4$.</figcaption></figure>

<table><caption>**Table 1.** $r^2$ values for the fit of $\ln(\epsilon)$ versus $\ln(\chi)$.</caption><thead><tr><th rowspan="2">q</th><th colspan="2">Simulated</th><th colspan="2">Real</th></tr><tr><th>φ = 0.400</th><th>φ = 0.08</th><th>φ = 0.399</th><th>φ = 0.073</th></tr></thead><tbody><tr><td>-9</td><td>0.986</td><td>0.986</td><td>0.898</td><td>0.935</td></tr><tr><td>-8</td><td>0.987</td><td>0.987</td><td>0.900</td><td>0.936</td></tr><tr><td>-7</td><td>0.989</td><td>0.989</td><td>0.902</td><td>0.938</td></tr><tr><td>-6</td><td>0.992</td><td>0.991</td><td>0.905</td><td>0.942</td></tr><tr><td>-5</td><td>0.994</td><td>0.993</td><td>0.910</td><td>0.946</td></tr><tr><td>-4</td><td>0.995</td><td>0.995</td><td>0.917</td><td>0.952</td></tr><tr><td>-3</td><td>0.996</td><td>0.997</td><td>0.928</td><td>0.962</td></tr><tr><td>-2</td><td>0.997</td><td>0.999</td><td>0.949</td><td>0.979</td></tr><tr><td>-1</td><td>0.998</td><td>0.999</td><td>0.984</td><td>0.998</td></tr><tr><td>0</td><td>0.998</td><td>0.999</td><td>0.996</td><td>0.999</td></tr><tr><td>2</td><td>0.998</td><td>0.999</td><td>0.995</td><td>0.999</td></tr><tr><td>3</td><td>0.998</td><td>0.999</td><td>0.995</td><td>0.999</td></tr><tr><td>4</td><td>0.999</td><td>0.998</td><td>0.995</td><td>0.999</td></tr><tr><td>5</td><td>0.999</td><td>0.998</td><td>0.995</td><td>0.999</td></tr><tr><td>6</td><td>0.999</td><td>0.998</td><td>0.996</td><td>0.999</td></tr><tr><td>7</td><td>0.999</td><td>0.998</td><td>0.996</td><td>0.999</td></tr><tr><td>8</td><td>0.999</td><td>0.998</td><td>0.996</td><td>0.999</td></tr><tr><td>9</td><td>0.999</td><td>0.998</td><td>0.996</td><td>0.999</td></tr></tbody></table>

high porosity. At high porosity, the curve remains quite linear, tending towards a monofractal structure. The deviation of the curve from linearity indicates that at low porosity, the sample shows multiple scaling behaviour, that is the low density and high density regions of pores scale differently.

Subsequently, the generalized dimension $D_q$ were estimated using eq. (12), and eq. (11) for $D_1$. Plots of $D_q$ versus $q$ for a low porosity corresponding to $p = 0.9$ and a high porosity value corresponding to $p = 0.8$ are displayed in Figs 8(a) and (b), respectively. For monofractals, all the $D_q$s would lie on the same horizontal line. The multifractal nature of the rock sample at low porosity is more pronounced than that of the higher porosity sample as evident from Figs 8(a) and (b). At high porosity, the sample has an

almost monofractal nature. The $D_q$ values corresponding to the first three moments and the interrelations between them for different porosities are displayed in Table 2. The difference in the values of the generalized dimension for the first three moments though small, is unmistakable, especially for low porosity, indicative of the multifractal nature of the samples.

The capacity dimension, $D_0$ provides information on how abundantly the measure, defined by eq. (3), is distributed over the scales of interest. From Table 2, the value of $D_0$ is less than the euclidean dimension 3 for the entire range of porosity studied. The $D_0$ is a maximum at the maximum porosity corresponding to $p = 0.5$ and has a minimum of 2.891 corresponding to $p = 0.9$.

The entropy or information dimension $D_1$ reaches a maximum of 2.925 corresponding to maximum porosity. Its value decreases on either side of the maximum porosity. Lower $D_1$ indicates greater concentration of pores over a small size domain, that is greater clustering. Comparison of Figs 2 and 3 clearly indicates that there is greater clustering of pores with decreasing porosity. At higher porosity, the pore clusters are more uniformly distributed which corresponds to a higher value of $D_1$. For distributions with a given value of $D_0$, the maximum possible value of $D_1$ is just $D_0$, which corresponds to the least heterogeneous distribution with such $D_0$. This suggests that the ratio $D_1/D_0$ may be seen as a measure of the dispersion of the porosity relative to the dispersion of the pore sizes. The closer to 1 this ratio is, the more evenly distributed is the porosity over the set of pore sizes.

The correlation function $D_2$ describes the uniformity of the measure (here pore cluster size). Smaller $D_2$ values indicate long-range dependence, whereas higher values indicate domination of short range dependence. From Table 2 we see that $D_2$ has a minimum at the lowest porosity which is indicative of long range correlations appearing between pores. The long range correlation is a manifestation of our growth algorithm for the sedimentary rocks. When the fraction of larger grains is small, elongated and isolated pore clusters are more prominent. Thus the pore sites show greater autocorrelation along these clusters.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Variation of mass exponent $\tau(q)$ versus $q$ for real and simulated rock structures. The open circles show data for simulated structure while ‘×’ indicate data of real sample. (a) Low porosity $\phi \simeq 0.07$. (b) High porosity $\phi \simeq 0.4$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Plot of $D - q$ versus $q$ for real and simulated rock structures. The open circles show data for simulated structure while ‘×’ indicate data of real sample. (a) Low porosity $\phi \simeq 0.07$. (b) High porosity $\phi \simeq 0.4$.</figcaption></figure>

<table><caption>**Table 2.** Comparison of generalized fractal dimensions for various probability values.</caption><thead><tr><th>p</th><th>φ</th><th>D(0)</th><th>D(1)</th><th>D(2)</th><th>D(1)/D(0)</th><th>D(0) - D(2)</th></tr></thead><tbody><tr><td>0.0</td><td>0.435</td><td>2.9190</td><td>2.9108</td><td>2.9042</td><td>0.99719</td><td>0.0148</td></tr><tr><td>0.1</td><td>0.443</td><td>2.9209</td><td>2.9126</td><td>2.9059</td><td>0.99716</td><td>0.0150</td></tr><tr><td>0.2</td><td>0.452</td><td>2.9226</td><td>2.9142</td><td>2.9074</td><td>0.99713</td><td>0.0152</td></tr><tr><td>0.3</td><td>0.458</td><td>2.9239</td><td>2.9154</td><td>2.9085</td><td>0.99709</td><td>0.0154</td></tr><tr><td>0.4</td><td>0.463</td><td>2.9248</td><td>2.9159</td><td>2.9089</td><td>0.99696</td><td>0.0159</td></tr><tr><td>0.5</td><td>0.464</td><td>2.9248</td><td>2.9157</td><td>2.9085</td><td>0.99689</td><td>0.0163</td></tr><tr><td>0.6</td><td>0.461</td><td>2.9236</td><td>2.9142</td><td>2.9068</td><td>0.99678</td><td>0.0168</td></tr><tr><td>0.7</td><td>0.449</td><td>2.9202</td><td>2.9103</td><td>2.9025</td><td>0.99661</td><td>0.0177</td></tr><tr><td>0.8</td><td>0.421</td><td>2.9123</td><td>2.9017</td><td>2.8932</td><td>0.99636</td><td>0.0191</td></tr><tr><td>0.9</td><td>0.356</td><td>2.8908</td><td>2.8789</td><td>2.8690</td><td>0.99588</td><td>0.0218</td></tr></tbody></table>

In an earlier study (Sadhukhan et al. 2007b) on conductivity through connected pore space of sedimentary rocks, the effective conductivity of the simulated rock using RBBDM, showed a maximum for $p = 0.7$ despite a maximum porosity at $p = 0.5$. For structures generated by using the RBBDM, the maximum backbone mass of the connected cluster corresponded to $p = 0.7$. The authors had established that it was the backbone mass of the connected cluster that was most effective for transport. The difference ($D_0 - D_2$) for different $p$ values are displayed in Table 2. As smaller $D_2$ values indicate better correlation, a higher value of ($D_0 - D_2$)

is indicative of a more homogeneous correlated pore distribution. Since a sedimentary rock should not only have sufficiently high porosity but also good correlation between the pores, we calculated the product $\xi = \phi(D_0 - D_2)$ for the different $p$ values. A variation of $\xi$ with $p$ is displayed in Fig. 9(a) along with the variation of conductivity $\sigma$ with $p$ (Sadhukhan et al. 2007b). The proposed structure factor follows a linear relationship with conductivity for different values of $p$ upto its maximum value, as displayed in Fig. 9(b). Beyond the maximum value, the plot curves backwards as expected for any extremum point. The linear relationship between $\xi$ and $\sigma$ leads us to propose that $\xi$ is the ‘structure parameter’ that determines conductivity in these porous rocks.

The $\alpha$ and $f(\alpha)$ values of the singularity spectrum were computed with the help of eq. (9). The Hölder exponent $\alpha$ for each $p$ gives the average values of local mass distribution for a given scale. The largest $f(\alpha)$ corresponding to $\alpha_0$, is obtained for $q = 0$. The maximum $f(\alpha)$ corresponds to the capacity dimension $D_0$. A greater value of $\alpha$ indicates a lower degree of mass concentration and vice versa. The parameters of the singularity curves shown in Fig. 10, are listed in Table 3 for the two different values of porosity. The convex hump of the singularity spectra for the two simulated rock samples clearly indicate the multifractal nature of the structure.

Width of the $f(\alpha)$ spectra is defined as the difference between the $\alpha$ values of the most negative moment $q_-$, that is $\alpha_{max}$, and the most positive moment $q_+$, that is $\alpha_{min}$. The wider the spectrum, that
<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) Plot of structure parameter $\xi$, and a plot of conductivity $\sigma$ with $p$ in simulated rocks using RBBDM. The two sets of data points show similar variation with respect to $p$. (b) Plot of structure parameter $\xi$ versus conductivity $\sigma$. The relationship is linear upto the maximum. The plot curves backwards beyond the maximum point.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** $f(\alpha)$ versus $\alpha$ for real and simulated rock structures. The open circles show data for simulated structure while ‘x’ indicate data of real sample. (a) Low porosity $\phi \simeq 0.07$. (b) High porosity $\phi \simeq 0.4$.</figcaption></figure>

<table><caption>**Table 3.** Comparison of local fractal dimensions.</caption><thead><tr><th rowspan="2">α</th><th colspan="2">Simulated</th><th colspan="2">Real</th></tr><tr><th>φ = 0.400</th><th>φ = 0.08</th><th>φ = 0.399</th><th>φ = 0.073</th></tr></thead><tbody><tr><td>α(0)</td><td>2.905</td><td>2.965</td><td>2.839</td><td>2.981</td></tr><tr><td>α<sub>max</sub></td><td>3.098</td><td>3.423</td><td>5.079</td><td>4.773</td></tr><tr><td>α<sub>min</sub></td><td>2.829</td><td>2.828</td><td>2.325</td><td>2.791</td></tr><tr><td>α<sub>max</sub> - α(0)</td><td>0.193</td><td>0.458</td><td>2.240</td><td>1.792</td></tr><tr><td>α(0) - α<sub>min</sub></td><td>0.076</td><td>0.137</td><td>0.514</td><td>0.190</td></tr><tr><td>α<sub>max</sub> - α<sub>min</sub></td><td>0.269</td><td>0.595</td><td>2.754</td><td>1.982</td></tr></tbody></table>

is greater the difference between ($\alpha_{max} - \alpha_{min}$), the higher is the heterogeneity in the scaling indices of pore mass and vice versa. As seen from Fig. 10(a), the width of the singularity spectrum is greater for the simulated rock structure of low porosity than the corresponding width of the spectrum of the high porosity simulated structure, inset of Fig. 10(b), in compliance with Figs 2 and 3. While Fig. 2 shows a homogeneity in pore size both in the bedding plane and its perpendicular plane at high porosity, the pore size and shape show greater heterogeneity in the corresponding planes for the low

porous sample. Small $f(\alpha)$ values indicate rare events (extreme values of the PSD). Asymmetry in the $f(\alpha)$ spectra indicate the dominance of higher or lower values of pore masses. A larger width of $\alpha_0 - \alpha_{min}$ indicates the domination of large values in the PSD. An examination of Table 3 shows that both $\alpha_0 - \alpha_{min}$ and $\alpha_{max} - \alpha_0$ decrease with decreasing $p$. On the other hand, a large right width $\alpha_{max} - \alpha_0$ would indicate the dominance of extremely small values in the PSD. The long tail of Fig. 5(a) for low porosity is reflected in the greater width of $\alpha_{max} - \alpha_0$, whereas the very narrow distribution of pore sizes of the more porous rock (Fig. 5 b), shows up as a very narrow width of ($\alpha_{max} - \alpha_{min}$). For the porosities studied (Table 3), for $p$ values ~0.9, $\alpha_{max} - \alpha_0$ becomes much bigger than $\alpha_0 - \alpha_{min}$. In this region clustering of pore sites into elongated isolated channels occur leaving larger sections of structure free of pores as also seen in Fig. 2.

## 3.2 Multifractal analysis of structure of real rocks

To compare our simulation results with real rock samples, X-ray tomography micrographs of 2-D sections of two real sedimentary rock
samples obtained from an oolitic limestone (pure calcite) from the Mondeville formation of Middle Jurassic age (Paris Basin, France), and a reefal carbonate from the Majorca Islands, Spain, were used. The limestone is composed of recrystallized oolite with a mean diameter of less than a few hundred µm. Each pixel of the micrographs 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. The binary file was then converted to a grey-scaled picture, as shown in Fig. 4, using ImageJ. 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 in size. Sections of the 3-D limestone rock cut along the direction of assembly are shown in Figs 4(a) and (b) while Fig. 4(c) shows a section of the same sample cut perpendicular to the direction of assembly. Similar sections were cut from the reefal carbonate and these are shown in Figs 4(d)–(f). The porosity of the oolitic limestone was determined from the reconstructed rock and found to be 0.074 while the reefal carbonate was found to have a porosity of 0.399. This high contrast in their porosity values is evident from the panels of Fig. 4. The PSD of the low and high porosity real rock samples were determined as before, and displayed in Figs 5(c) and (d), respectively. The low porous limestone clearly is more heterogeneous in size distribution than the more porous carbonate rock. However both the real samples have a long tail distribution.

The plots of $\chi(q, \epsilon)$ versus $\epsilon$ for both the real rock samples are displayed in Figs 6(c) and (d), while the coefficient of determination ($R^2$) for the same samples are displayed in Table 1. The $R^2$ values for real rock of high porosity become less than 0.91 as $q < -5$. However to enable a comparison between real and simulated rocks, we have carried out our analysis over the same range of $q$ values as the simulated rock. The variation of $\tau(q)$, $D_q$ versus $q$ and the $f(\alpha)$ spectra for the real rocks have been plotted along with their closest matching porosity samples generated by the RBBDM.

The non-linear nature of the plots $\tau(q)$ versus $q$, with two distinct slopes for positive and negative $q$ values clearly indicate that both the real rock samples are also multifractal (Figs 7a and b). The match between real and simulated sample for low porosity is remarkable except for negative $q$ values. Here of course the moment calculation emphasizes the regions of less pores. The high porosity carbonate sample has a more pronounced non-linear nature indicating an emphasized multifractal character.

The real rocks show a significant variation of $D_q$ with respect to $q$, having a nearly sigma shape (Fig. 8). The quasi-linear shape for positive $q$ values may be interpreted as a rather homogeneous distribution of high concentrations of porosity along the interval of pore sizes. The convex shape of the curves for negative $q$ values indicate a heterogeneous distribution of low concentrations of porosity. Once again the match between the real and the simulated structure for low porosity and positive moments, is almost exact. The mismatch between real and simulated rocks for high porosity show that though the RBBDM can generate 3-D porous structures of required porosity, the PSD and correlation between the pore clusters is not captured at high porosities. Clearly the real rocks are more multifractal with the $D_q$ values decreasing rapidly with higher moments.

The multifractal spectrum of the limestone rock indicates greater dominance of larger pores and shows less heterogeneity in the $f(\alpha)$ values (Fig. 10a). The more porous carbonate sample shows a symmetric distribution of $f(\alpha)$ values with greater heterogeneity of $f(\alpha)$ for larger clusters (Fig. 10b). The match in the $f(\alpha)$ spectrum be-

tween the real and simulated structures for low porosities is almost exact.

## 4 CONCLUSIONS

The PSD and microgeometry of the pore phase of a sedimentary rock can be quite different when examined at different length scales. The length scale of interest depends on the process being studied. Multifractal analysis not only shows up the scaling of a measure at different length scales, but the generalized dimensions give additional information of the average, the concentration or clustering of the measure and correlation between the measures as we calculate higher moments. For a two-phase system, apart from information about the distribution of the measure of interest, the analysis gives us information on the complement of the measure at the same time.

We have done a multifractal analysis on sedimentary rock structures simulated by using RBBDM. The structures at different porosities all showed multifractal characteristics but the heterogeneity of the multifractal spectrum decreases with higher porosities. To have an idea of the suitability of RBBDM as a realistic model for generating sedimentary rocks, we carried out similar analysis on two real sedimentary rocks which have almost identical porosity values. The low porous limestone rock shows an almost identical match in its multifractal parameters with its corresponding simulated structure. The high porosity carbonate sample was however, quite different from the structure generated by RBBDM with identical porosity. Visually, sections of real and simulated samples seemed to indicate that the similarity was greater between the two at high porosity values. Comparison through multifractal analysis proved the contrary—the simulated samples actually have an almost match of their generalized dimensions at lower porosities than at higher values.

Transport properties of sedimentary rocks are expected to depend on the porosity as well as PSD. It is expected that the generalized dimensions of a rock should give information on transport properties of the rock as well. This study has lead us to propose a structure parameter $\xi$ which is a function of porosity and the generalized dimensions of multifractal rock, which effectively controls conductivity of a rock structure. Further, the structure parameter follows a linear relation with respect to the conductivity of the rock sample till conductivity is optimized. The optimal conductivity, guided by the optimization of the structure factor, may depend on the size and shape of the grains. Further studies with grains of different aspect ratio and sizes are required in order to understand the interdependence of optimal conductivity and structure factor. We also plan to test the validity of our proposition by studying other transport properties of simulated and real rocks in context with varying porosity and generalized dimensions.

## ACKNOWLEDGEMENTS

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

Chhabra, A.B., Meneveau, C., Jensen, R.V. & Sreenivassen, K.R., 1989. Direct determination of the $f(\alpha)$ singularity spectrum and its application to fully developed turbulence, *Phys. Rev. A*, **40**, 5284–5294.
Dathe, A., Tarquis, A.M. & Perrier, E., 2006. Multifractal analysis of pore and solid phases in binary two-dimensional images of natural porous structures, *Geoderma*, **134**, 318–326.

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

Feder, J., 1988. *Fractals*, Plenum Press.

Giri, A., Tarafdar, S., Gouze, P. & Dutta, T., 2012. Fractal pore structure of sedimentary rocks: simulation in 2-D using a relaxed bidisperse ballistic deposition model, *J. appl. Geophys.*, **87**, 40–45.

Grassberger, P. & Procaccia, I., 1983. Characterization of strange attractors, *Phys. Rev. Lett.*, **50**(5), 346–349.

Grau, J., Mendez, V., Tarquis, A.M., Diaz, M.C. & Saa, A., 2006. Comparison of gliding box and box-counting methods in soil image analysis, *Geoderma*, **134**, 349–359.

Halsey, T.C., Jensen, M.H., Kadanoff, L.P., Procaccia, I. & Shraiman, B.I., 1986. Fractal measures and their singularities: the characterization of strange sets, *Phys. Rev. A*, **33**(2), 1141–1151.

Hansen, J.P. & Skjeltorp, A.T., 1988. Fractal pore space and rock permeability implications, *Phys. Rev. B*, **38**, 2635.

Hansen, J.P., McCauley, J.L., Muller, J. & Skjeltorp, A.T., 1988. Multifractal analysis of sedimentary rocks—random fluctuations and pattern growth: experiments and models, *NATO ASI Series*, **157**, 310.

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.

Lovoll, G., Jankov, M., Mly, K.J., Toussaint, R., Shmittbuhl, J., Schfer, G. & Mheust, Y., 2011. Influence of viscous fingering on dynamic saturation-pressure curves in porous media, *Transport Porous Media*, **86**(1), 305–324.

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

Nsir, K., Schfer, G., Roupert, R.D.C., Razakarisoa, O. & Toussaint, R., 2012. Laboratory experiments on DNAPL gravity fingering in water-saturated porous media, *Int. J. Multiphase Flow*, **40**, 83–92.

Paz-Ferreiro, J., Wilson, M. & Vidal Vázquez, E., 2009. Multifractal description of nitrogen adsorption isotherms, *Vadose Zone J.*, **8**, 209–219.

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

Rényi, A., 1960. On measures of entropy and information, in *Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability*, Vol. 1 (University of California Press, Berkeley, 1961), pp. 547–561.

Rieu, M. & Sposito, G., 1991. Fractal fragmentation, soil porosity and soil water properties, 1. Theory, *Soil Sci. Soc. Am. J.*, **67**, 1361–1369.

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

Shannon, C., 1949. Communication Theory of Secrecy Systems, *Bell Sys. Tech. J.*, **28**(4), 656–715.

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

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

Tarquis, A.M., Gimenez, D., Saa, A., Diaz, M.C. & Gasco, J.M., 2003. Scaling and multiscaling of soil pore systems determined by image analysis, in *Scaling Methods in Soil Physics*, pp. 19–34, eds Selim, H.M., Pachepsky, Y. & Radcliffe, D.E., CRC Press.

Tarquis, A.M., Heck, R.J., Grau, S.B., Fabregat, J., Sanchez, M.B. & Anton, J.M., 2007. Influence of thresholding in mass and entropy dimension of 3-D soil images, *Nonlin. Process. Geophys.*, **15**, 881–891.

Toussaint, R., Lovoll, G., Mheust, Y., Mly, K.J. & Schmittbuhl, J., 2005. Influence of pore-scale disorder on viscous fingering during drainage, *Europhys. Lett.*, **71**, 583–589.

Vicsek, T., 1992. *Fractal Growth Phenomena*, 2nd ed., World Scientific.

Vidal Vázquez, E., Paz-Ferreiro, J., Miranda, J.G.V. & Gonzalez, A.P., 2008. Multifractal analysis of pore size distributions as affected by simulated rainfall, *Vadose Zone J.*, **7**, 500–511.