Estimation of permeability of a sandstone reservoir by a fractal and Monte Carlo simulation approach : a case study

Permeability of a hydrocarbon reservoir is usually estimated from core samples in the laboratory or from well test data provided by the industry. However, such data is very sparse and as such it takes longer to generate that. Thus, estimation of permeability directly from available porosity logs could be an alternative and far easier approach. In this paper, a method of permeability estimation is proposed for a sandstone reservoir, which considers fractal behavior of pore size distribution and tortuosity of capillary pathways to perform Monte Carlo simulations. In this method, we consider a reservoir to be a mono-dispersed medium to avoid effects of micro-porosity. The method is applied to porosity logs obtained from Ankleshwar oil field, situated in the Cambay basin, India, to calculate permeability distribution in a well. Computed permeability values are in good agreement with the observed permeability obtained from well test data. We also studied variation of permeability with different parameters such as tortuosity fractal dimension ( Dt), grain size ( r) and minimum particle size ( d0), and found that permeability is highly dependent upon the grain size. This method will be extremely useful for permeability estimation, if the average grain size of the reservoir rock is known.


Introduction
Permeability is one of the important parameters that govern production of hydrocarbons from reservoirs.As we know, fluid flow in reservoirs depends upon permeability, and in the case of multiple phases of hydrocarbons, it is relative permeability that governs the fluid flow.Estimation of the permeability of a geological medium is a difficult job, as it varies over several orders of magnitude, even for a single rock such as sandstone (Clauser, 1992;Nelson, 1994;Dimri et al., 2012).In the literature many empirical relations are available for computation of permeability from porosity (Kozeny, 1927;Carman, 1956;Pape et al., 1999).These empirical relations define the porosity-permeability relationship in terms of correlation coefficients without addressing the real physical situation existing in a porous medium.The newly developed nuclear magnetic resonance tool can be used to estimate reliable permeability, but it is very new and still not included in the commonly recorded log suite.In general, permeability is estimated from effective porosity.However, even for a given effective porosity, permeability will be different for different rock types.This is because in addition to porosity, permeability depends upon textural parameters such as grain size, grain shape and sorting of grains.These textural parameters affect the tortuous nature of capillary pathways and the arrangement of pores in the porous medium, thereby making it more complex to understand.The concept of fractals is useful for solving many complex problems in the earth sciences (Dimri, 2005;Dimri et al., 2005;Dimri et al., 2012;Vedanti and Dimri, 2003;Vedanti et al., 2011;Srivastava and Sen, 2010, etc.) and fractals have also been useful for explaining the complex nature of porous media (Adler and Thovert, 1993;Dimri, 2000a, b;Feranie and Latief, 2013;Katz and Thompson, 1985;Krohn and Thompson, 1986;Pape et al., 1999;Smidt and Monro, 1998;Young and Crawford, 1991;Yu and Li, 2001, etc.).There are numerous examples that present the application of fractal geometry to analyze porous media.Sahimi and Yortsos (1990) presented a review of the general classes of application of fractal geometry to porous media.Without involving any empirical constants, Pitchumani and Ramakrishnan (1999) expressed permeability in terms of two fractal dimensions: one describes the size distribution of capillary pathways and the other describes the tortuosity of the capillary pathways.2001) deduced a unified model for describing the fractal character of porous media and proposed a criterion to determine whether a porous medium can be characterized by fractal theory and technique or not.An analytical expression for the permeability of a fractal bi-dispersed porous medium was developed by Yu and Cheng (2002).However, if a porous medium is a multiscale or scale-dependent fractal medium, analytical expression may not give accurate results.To overcome this issue, Yu et al. (2005) proposed the Monte Carlo technique to simulate the permeability of fractal porous media.Yu and Li (2004) derived analytical expressions of the fractal dimensions for wetting and non-wetting phases for unsaturated porous media.Liu and Yu (2007) established a fractal relative permeability model that takes into account the capillary pressure difference effect in the case of unsaturated porous media.Yu (2008) reviewed the theories and achievements in the field of application of fractal geometry to understand flow in porous media.Xu and Yu (2008) derived an expression for the Kozeny-Carman constant by using the analytical formula of permeability.The Monte Carlo simulation technique is used by Xu et al. (2012) to estimate relative permeability.A fractal model for capillary pressure is obtained by Xiao and Chen (2013).Xu et al. (2013a) used fractal theory and the Monte Carlo simulation technique to develop a probability model for radial flow in fractured porous media.Analytical expressions for relative permeabilities of wetting and non-wetting phases are presented by Xu et al. (2013b).The innovative work presented by Yu et al. (2005) assumed a model in which a porous medium consists of a set/bundle of parallel and tortuous capillaries with uniform diameter and obtained probability model for pore diameter and permeability.The probability model comprehends the fractal nature of pore size distribution and tortuous capillaries.The overall approach consists of simulating random pore size distribution to calculate flow rate from the Hagen-Poiseuille equation (Denn, 1980) and then, by using Darcy's law (Darcy, 1856), an expression for permeability is obtained.This fractal-based method was applied to estimate the permeability of a sintered copper bi-dispersed porous medium to solve a heat transport problem (Yu et al., 2005).
The estimation of the permeability of a hydrocarbon reservoir is one of the important goals of reservoir geophysicists.In the present work, a fractal-based Monte Carlo simulation approach is applied for the first time to estimate the permeability of the Ankleshwar sandstone hydrocarbon reservoir, situated in the Cambay basin, India.The reservoir formation consists mainly of alternate layers of sandstone and shale.As we know that sandstones and shales are porous fractal media (Katz and Thomson, 1985;Krohn and Thompson, 1986;Krohn, 1988b;Sahimi, 2011), we performed Monte Carlo simulations based on the fractal nature of pore size distribution in porous media of the reservoir.We assumed the reservoir as a mono-dispersed porous medium to study the effect of macro-porosity on permeability.A mono-dispersed porous medium considers only macro-pores formed between the grains and ignores micro-pores inside the clusters, which do not contribute towards macro-porosity.

Fractal description of pore structures
A porous medium consists of tortuous capillary pathways.The tortuous length of these capillary pathways follows the fractal/scaling law.The scaling relationship for capillary length in heterogeneous porous media is given by L t (ε) = ε 1−D t L D t 0 , where ε, L t and L 0 are the scale of measurement, tortuous length and characteristic (straight) length, respectively and D t is the tortuosity fractal dimension (Wheatcraft and Tyler, 1988).Yu and Cheng (2002) considered pore diameter (λ) as the scale of measurement.Thus, the smaller the diameter, the longer the capillary and scaling law that can be given as where 1 < D t < 2, representing the extent of convolutedness of capillary pathways for fluid flow through a medium.The higher the value of D t , the more the convolutedness or tortuosity.The limiting values of D t = 1 and D t = 2 correspond respectively to a straight capillary and a highly tortuous line that fills a plane (Wheatcraft and Tyler, 1988).
A porous medium consists of a large number of pores of varying pore diameters that intersect the pore cross sections.The size distribution of pore diameters is another important property.This distribution is analogous to cumulative size distribution of islands on the Earth's surface, which follows the fractal scaling law (Mandelbrot, 1982;Majumdar and Bhushan, 1990).Pitchumani and Ramakrishnan (1999) and Yu and Cheng (2002) have established a fractal scaling law to describe the distribution of pores in a porous medium, which is given as: where N is the number of pores with diameter (L) greater than or equal to λ, D f the pore area dimension, with 1 < D f < 2, representing the fractal dimension of the intersecting pore cross sections with a plane normal to the flow direction.It is evident from Eq. ( 2) that when λ approaches maximum pore size, λ max , the number of pores greater than or equal to λ max is one.Conversely, when λ approaches smallest pore size, λ min , the numbers of pores are maximum, and the scaling law becomes where N t is the total number of pores.
Based on these fractal scaling laws, Yu et al. (2005) obtained probability expressions for pore diameter and permeability in terms of pore diameter as given below: where (7) In the case of a bi-dispersed porous medium, effective porosity is given by Yu and Cheng (2002) as below.
where φ e = effective porosity, φ c = micro-porosity inside cluster, R c = cluster mean radius, d + = ratio of cluster mean diameter to the minimum particle size (d 0 ), λ i = diameter of the ith capillary tube chosen by the Monte Carlo simulation, A = total cross-sectional area of a unit cell, and G = geometry factor for flow through a circular capillary.
In a bi-dispersed porous medium, small particles form clusters.The spaces between clusters form macro-pores and within the cluster micro-pores exist.

Geology of Ankleshwar oil field, Cambay basin
Ankleshwar oil field lies in the western part of the Cambay basin, India (Fig. 1).The field is doubly plunging anticline, trending ENE-WSW.It has a multi-layered sandstone reservoir of deltaic origin.Ankleshwar formation, which was deposited during the marine regression phase during the Upper to Middle Eocene (Holloway et al., 2007) is the major stratigraphic unit in the field.The formation consists of four sub-lithological units, viz., Telwa shale (effective seal), Ardol, Kanwa shale and Hazad members from top to bottom (Fig. 5).A Hazad member is an important reservoir sandstone; however, it contains shale laminae.The Ardol section comprises sand-shale alterations.The Ankleshwar formation is overlain by the Dadhar formation and underlain by Cambay shale, which is the source rock.

Estimation of porosity from density log
Porosity is derived from density log using the following equation: where ρ b = bulk density of the formation, ρ m = density of the rock matrix, ρ f = density of the fluids occupying the porosity, and φ = porosity of the rock.In Eq. ( 9) the value of ρ m is taken as 2.65 g cm −3 , which is the default density of quartz grains, and the value of ρ f is taken as 1.1 g cm −3 , which is the default density of saline water.Equation ( 9) corresponds to porosity of pure sandstone.Thus, to discriminate between pure sand, pure shale and shaly-sand intervals in the well, the criteria shown in Table 1 were adopted.We know that shale is more porous than sandstone, but pores are not interconnected, and hence cannot permeate fluid.Hence, to find effective or interconnected porosity of shaly-sand intervals, bulk density (ρ b ) is corrected for shale volume using Eq. ( 10), which is given as: where (ρ b ) corr = corrected bulk density of the shaly lithology, (ρ b ) clean = bulk density of clean sand stone formation, (ρ b ) shale = bulk density of pure shale, and V sh = shale volume, which is calculated from the gamma-ray log using the Dresser tertiary equation.The values of (ρ b ) clean and (ρ b ) shale are 2.19 g cm −3 and 2.36 g cm −3 , respectively, which are selected from depth intervals 1150 m to 1158 m and 1184 m to 1194 m, respectively, in the log.The corrected value of density is then used in Eq. ( 9) to obtain effective porosity of shaly-sand intervals, which is given as: In shaly-sand intervals clay minerals occupy pores or they cover the sand grains, thereby forming micro-pores and introducing micro-porosity (φ c ), which is the porosity below which there is no permeating flow (the percolation threshold) (Nabovati et al., 2009).In such cases, in place of effective porosity, useful porosity (φ use ) (http://www.spec2000.net)can be used, which is given as: where C is a constant that can be obtained from shale volume, as shown below: Here, we assume that permeability estimation by using useful porosity is a more efficient way, as in this case the contribution of micro-porosity to permeability becomes negligible.In the present work, a well W1 drilled in the Ankleshwar sandstone reservoir cutting across the Telwa, Ardol, Kanwa, Hazad and Cambay lithological units is used for the analysis.The Hazad and Ardol sections of the reservoir are divided into eight sand layers, which are named S1 to S8. Intervals of these sand layers with useful porosity greater than 15 % are considered as reservoir intervals.The average useful porosity obtained from the reservoir intervals for each sand layer is given in Table 3.The calculated useful porosity of different sand layers matches with the information given by the operator in this well.

Monte Carlo simulations for prediction of permeability
Monte Carlo simulation is a random search method, widely used in geophysics, using generation of random numbers (Dimri, 1992).The method allows the generation of many possible models.In the present work, the method suggested by Yu et al. (2005) is applied to the Ankleshwar sandstone reservoir, where micro-pores are formed by pore-filling material (clays) or exist as intra-granular pores (Krohn, 1988a).
The permeability contribution by these micro-pores is negligible (Yu and Lee, 2000;Nimmo, 2004;Loucks, 2005).Hence, as discussed above, it is more appropriate to run the simulations with macro-porosity or φ use .After removal of micro-porosity, clusters become solid grains; hence instead of cluster mean radius, grain radius (r) can be used.In this case Eq. ( 8) becomes: The flow chart to implement the method is shown in Fig. 2 and the algorithm for determination of the permeability from macro-porosity is given in Appendix A.

Variation of permeability with tortuosity fractal dimension (D t ), and average grain radius (r) and minimum particle size (d 0 )
Monte Carlo simulations are run for different values of grain radius, tortuosity fractal dimension and minimum particle size to understand permeability variation with these parameters.The mean value of increment in permeability for a small variation in grain radius from 0.05 mm to 0.1 mm is 38 %, as shown in Fig. 3a (when D t = 1.8 and d 0 = 1 µm).Similarly, when D t = 1.8 and r = 0.05 mm, the minimum particle size (d 0 ) is increased from 1 µm to 2 µm, then permeability is increased by 22 %, as shown in Fig. 3b.Next, when D t is increased from 1.6 to 1.8 by considering d 0 = 1 µm and r = 0.05 mm, permeability is reduced by 20 %, as shown in Fig. 3c.In each case D f is calculated from Eq. ( A4).The estimated permeability for different values of r, D t , d 0 and D f is given in Table 2.This analysis suggests that permeability is more sensitive to change in grain radius than D t and d 0 .

Estimation of grain radius
Sand grains in different reservoir layers can have different radii.Mavko and Nur (1997) incorporated micro-porosity into the Kozeny-Carman equation, which we used to compute the distinct grain radius of reservoir layers and the equation is given below: Equation ( 17) is based on the assumption that formation is made up of spherical grains of uniform diameter.
However, we have assumed that the observed average reservoir permeability and average porosity pairs as shown in Table 3, available from well test data, are those of pure sandstone layers because the reservoir consists of a negligible amount of shale.In this case φ c is zero, φ e is equal to φ use and thus in this case grain radius, r will become r sand , which is the radius of sand grains and can be given as where K resvr is the observed average reservoir permeability.In Eq. ( 18) tortuosity (τ ) is unknown.The average grain radius of the S 1 and S 3 layers measured from core samples in nearby wells is 0.069 mm and 0.142 mm, respectively (R. Sharma, personal communication, 2012).Thus, using known grain radii, the average value of τ is estimated as 6.99.For other sand layers (S 2 , S 4 to S 8 ), r sand is estimated from Eq. ( 18) and the values are given in Table 3.However, in case of shaly-sand intervals, the effective grain radius (r eff ) is calculated using a weighted average formula using volume fractions as the weights, which is given as

Assigning the values to tortuosity fractal dimension (D t ) and minimum particle size (d 0 )
In order to select the values of D t and d 0 , simulations are run for several values (1.75, 1.8, 1.85 and 1.87) of D t and three different values (1.0 µm, 1.5 µm, 2.0 µm) of d 0 , with φ use (Fig. 5) and D f (Eq.A4) as the known inputs.Since litholog given by the operator shows the presence of clay in all the formations, we assumed minimum particle size (d 0 ) of clay grade, which are defined as the grains with 1 to 2 µm diameters.As we know that permeability decreases with an increase in D t (Fig. 3c) and observed reservoir permeability values are low, we chose higher values of D t for the analysis.For all possible pairs of D t and d 0, , permeability of all the sand layers is calculated.In each sand layer the average of useful porosity and permeability over all the reservoir intervals (Table 4) is calculated and compared with the corresponding observed average porosity and average reservoir permeability available from well tests (Table 3).Root mean square error (RMSE) between calculated and observed values is measured for each pair of D t and d 0 .

Results and discussions
The calculated average reservoir permeability of each sand layer for all the pairs of D t and d 0 and the corresponding RMS error are given in Table 4.It is clear from Table 4 that D t = 1.85 and d 0 = 1 µm give the least value of the RMS error, which is 53.4 mD.The mean value of D f for all sand layers, calculated from Eq. (A4), is 1.67.The RMS error between observed and calculated porosity is 1.8 %.Permeability of all sand layers is also calculated from the Kozeny-Carman equation (Eq.17) by replacing φ e − φ c with φ use and "r" with r eff .In this case the RMS error between calculated and observed values is 119.6 mD (Table 4).
The values of tortuosity fractal dimension, minimum particle size and pore area fractal dimension used to calculate the permeability log in the Ankleshwar formation are 1.85, 1 µm and 1.67, respectively.The calculated reservoir permeability versus observed reservoir permeability plot is shown in Fig. 4. The blocked logs of volume of sand, volume of U. Vadapalli et al.: Estimation of permeability of a sandstone reservoir shale, effective porosity, critical porosity, useful porosity and permeability are shown in Fig. 5.The highest value of permeability (409 mD) is obtained for reservoir zone S3, which is highlighted by rectangles in Fig. 5.This value exactly matches with the observed permeability of the S3 layer obtained from well test data.For the S4 sand layer, the calculated value of permeability (27 mD) is very close to the observed non-reservoir permeability (23 mD) provided by the operator.
The permeability of each sand layer in the Ankleshwar formation is estimated with the confirmed set of fractal dimensions, minimum particle size and estimated grain radius.Figure 4 clearly shows that, for most of the sand layers, the error between calculated and observed permeability does not exceed 50 mD, except in the case of the S1 sand layer, where the error is 86 mD.As we know, permeability obtained from well test data corresponds to a larger area in reservoir; however, permeability obtained by logs represents a smaller area around the well.This may also result in a mismatch between observed and calculated permeabilities.Another reason for the difference between observed and calculated permeabilities can be grain size, which is a very important parameter.
Since most of the evaluated permeability values match with the corresponding observed values within acceptable error range, the estimated values of D t and d 0 are reliable.From Fig. 5 it is obvious that porosity and permeability logs clearly discriminate pure sand, shaly sand and shale intervals, with higher values for pure sand, lower values for shaly sand and lowest values for shale.In the present well, S4 sand having very low permeability (27 mD) is non-reservoir, which is because of silt stone present in this interval (according to the litholog given by the operator).The Kanwa interval is defined as shale according to general geology; however, in the present well, according to the litholog given by the operator, a small amount of sandstone present in this interval reduces gamma-ray readings and thus this interval is shaly sand.It is obvious from the RMS error values given in Table 4 that Monte Carlo simulations give better results than the Kozeny-Carman equation.
Thus, by considering a mono-dispersed porous medium and modifying the method given by Yu et al. (2005), we could develop a methodology for estimation of reservoir permeability without using any empirical constant.The flow chart of the modified algorithm is also presented.

Conclusions
We found that accurate permeability estimation requires a good control of grain radius.In the major pay zone of the reservoir, the calculated permeability value (409 mD) exactly matches with the observed permeability value provided by the operator.In other sand layers, the difference between calculated and observed average reservoir permeability is ±50 mD, which is acceptable in this case.
Thus, in the absence of well test data or laboratory measurements, this method can be used to obtain first-hand information on reservoir permeability.Using this method, we can obtain continuous permeability distribution in reservoirs if porosity distribution from seismic data is known.9. Whenever convergence criterion is satisfied, calculate the permeability K (with D t measured from the box-counting method or using the method given in Sect.3.5).
10. Repeat the procedure for N number of times (say N = 1000) to get N values of permeability and take the mean to compute the final value of permeability.
11. Check the criterion, K min < K < K max , if not satisfied return to step 4 (where, K min and K max are acceptable minimum and maximum reservoir permeabilities respectively).

Fig. 1 .
Fig. 1.Location of the Cambay basin and its oil and gas fields (after www.spgindia.org/paper/sopt_2313/tmp_2313).Ankleshwar field is highlighted by the ellipse in the figure.

Fig. 2 .
Fig. 2. Flow chart to apply the Monte Carlo simulation technique for permeability estimation.
Fig. 3. (a), (b) and (c), respectively shows permeability variations grain radius (r), minimum particle size (d 0 ) and tortuosity fractal dimension (D t ).Permeability decreases with an increase in D t and increases with an increase in r and d 0 .The variation in permeability is more sensitive to changes in grain radius.

Fig. 4 .
Fig. 4. Permeability versus porosity plot for well (W1).Calculated average reservoir permeability (circles) falls in the range of observed average reservoir permeability values (stars).Sand layer names (S1 to S8) annotated to the data points clearly show the amount of deviation from predicted permeability with respect to the corresponding observed value.

Fig. 5 .
Fig. 5. Well logs of well (W1).Pure sandstone (< 35 API) and shale (> 65 API) are colored in red and blue on the gamma ray curve, respectively.Pink color bars and green color bars represent pure sandstone and shaly sand (≥ 35 API or ≤ 65 API), respectively.The zone highlighted by rectangles corresponds to a pure sandstone interval in the S3 layer that shows the highest permeability (409 mD), which matches with the observed value.

Yu Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union. U. Vadapalli et al.: Estimation of permeability of a sandstone reservoir and
Li (

Table 1 .
Gamma ray values adopted to discriminate pure sand, pure shale and shaly sand.

Table 2 .
Estimated permeability (K) for a range of porosities and for different values of grain radius (r), minimum particle size (d 0 ), tortuosity fractal dimension (D t ) and pore area fractal dimension (D f ).

Table 3 .
The observed average useful porosity, observed average reservoir permeability (k resvr ), calculated average useful porosity of reservoir intervals and estimated grain radius for different sand layers in the formation.

Table 4 .
Average values of reservoir permeability in each layer estimated for different values of D t and d 0 and from the Kozeny-Carman equation.For each estimated set the RMS error with respect to observed average reservoir permeability is given.