Distribution of Petrophysical Properties for Sandy-clayey Reservoirs by Fractal Interpolation

The sandy-clayey hydrocarbon reservoirs of the Upper Paleocene and Lower Eocene located to the north of Veracruz State, Mexico, present highly complex geological and petrophysical characteristics. These reservoirs, which consist of sandstone and shale bodies within a depth interval ranging from 500 to 2000 m, were characterized statistically by means of fractal modeling and geostatistical tools. For 14 wells within an area of study of approximately 6 km 2 , various geophysical well logs were initially edited and further analyzed to establish a correlation between logs and core data. The fractal modeling based on the R/S (rescaled range) methodology and the interpolation method by successive random additions were used to generate pseudo-well logs between observed wells. The application of geostatistical tools, sequential Gaussian simulation and exponential model var-iograms contributed to estimate the spatial distribution of petrophysical properties such as effective porosity (PHIE), permeability (K) and shale volume (VSH). From the analysis and correlation of the information generated in the present study, it can be said, from a general point of view, that the results not only are correlated with already reported information but also provide significant characterization elements that would be hardly obtained by means of conventional techniques .


Introduction
In the last years (Mandelbrot, 1983;Korvin, 1992;Barton and La Pointe, 1995), fractal geometry and its concepts have been considered as essential tools in many areas of the natural sciences (Vicsek et al., 1994;Sornette, 2006) due to the fact that the variation of the properties of many physical systems displays a fractal character.The thickness of lacustrine sediments, geological sediment logs and annual flood cycles in most rivers, for example, have exhibited long interdependence periods (Daryin and Saarinen, 2006).Hence, it is reasonable to expect a fractal character in the distribution of sediments because their statistics is determined by the natural processes that created them.In recent years, the concepts regarding fractal geometry have been applied for modeling the heterogeneity of reservoirs (Hewett, 2001;Srivastava and Sen, 2009).Applying fractal geometry to the description and assessment of reservoirs has a solid basis (Vivas, 1992;Yeten and Gümrah, 2000;Srivastava and Sen, 2010) since the distribution of the properties in sedimentary environments shows a fractal behavior with long-range correlations.Thus, fractal simulations can be used to generate distribution models of petrophysical and geological properties.
In the present research, sandy-clayey hydrocarbon reservoirs located to the north of Veracruz State, Mexico, with highly complex geological and petrophysical characteristics, and consisting of alternate sandstone and shale bodies (lutites) from the Upper Paleocene and Lower Eocene within a depth interval ranging from 500 to 2000 m, were studied by means of fractal modeling and geostatistical tools.The study area extends 123 km in length and 25 km in width (Abbaszadeh et al., 2003;Takahashi et al., 2006) and displays a set of submarine fans and turbiditic sediment deposits in a deep-water eroded canyon.These sediments show important variations concerning their clay-shale content in addition to altered secondary porosity due to diagenesis (Bermúdez et al., 2006;Talwani, 2011).The most important challenges that these reservoirs represent with respect to the improvement of their statistical characterization are the modeling of Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.petrophysical properties such as effective porosity (PHIE), which is the pore space from which fluids can be produced, permeability (K) and shale volume (VSH) at different well locations.The average neutron porosity and density were obtained and corrected by shale volume.
In order to fulfill this goal, based on both gamma-rayand porosity logs and well cores, a suitable lithologic model has been obtained for the study of the fractal characteristics and the determination of the petrophysical properties.The analysis of processes and interpretation of geophysical logs such as caliper (CALI), spontaneous potential (SP), gamma ray (GR), resistivity (LLD, ILD, MSFL), density (RHOB), neutron porosity (NPHI), sonic, and in some cases, water saturation (Sw) and permeability (K) were carried out for 14 wells.As for the improvement of the modeling of properties between wells, where data is not available, cross sections based on pseudo-well logs were obtained through fractal interpolation between neighboring well logs.The Hurst coefficient, which is necessary to perform this interpolation, was obtained by means of the rescaled range method (Hurst, 1951;Hurst et al., 1965) and applied to the geophysical well logs.In the present work, the compilation and analysis of data are presented, including the geological model (Talwani, 2011), along with the interpretation of petrophysical data used during the fractal characterization of the reservoirs (Arizabalo et al., 2004;Oleschko et al., 2008).From the analysis and correlation of the information generated in the present study, it has been found that the results not only are correlated with already reported information but also provide significant characterization elements that would be hardly obtained by conventional techniques.The local predictions regarding the high porosity and permeability, and low shale volume, represent a relatively high concentration of hydrocarbons in the area of study.

Statistical characterization of sandy-clayey reservoirs applying fractal and geostatistics modeling
By applying the R/S rescaled range methodology (Feder, 1988;Korvin, 1992;Srivastava and Sen, 2009) and the method of successive random additions (Voss, 1985(Voss, , 1988;;Saupe, 1988), a lateral interpolation was carried out to generate pseudo-well logs between observed well logs, showing, in addition, the necessary steps to perform the statistical characterization of a reservoir by fractal modeling.The application of geostatistics and fractal geometry involves the following steps: Selection of the reservoir, geophysical well logs and reservoir cores; assessment of the geological frame and complementary geophysical information; location and possible well connections; typical reservoir variogram (spatial variability function); fractal interpolation or well stochastic studies; pseudo-well petrophysical solution; vertical variation of well properties; cross sections of the porosity and permeability variations; identification and distribution of the reservoir flow units; representative variograms (flow units); areal distribution of the petrophysical parameters (flow units); and charts of the variation of the reservoir petrophysical parameters (Vivas, 1992).
In order to apply the methodology described above, it was necessary to verify the fractal behavior of the well logs; in order to do so, a test of the characteristics concerning the fractional Gaussian noise of the well logs was considered.Figure 1 shows the shale volume traces for the 14 studied wells, where superficial zones with high shale contents and deeper zones, where the shale volume is lower and the presence of hydrocarbons has been detected, can be seen.The VSH2 trace was chosen, which corresponds to the shale volume of well 2, to perform the fractal behavior test (Fig. 2).By means of the software BENOIT TM , the aforementioned test was carried out, which reflected an fGn type behavior (fractional Gaussian noise).The Rescaled Range method (R/S), Power Spectrum, Roughness-Length and Variogram were applied to the normalized trace values, showing, all of them,   a power law behavior (fractal).This procedure was applied to the 14 wells in the area of study.By comparing the H values obtained with the Rescaled Range and Roughness-Length methods, it was found that they are in good agreement, 0.5<H<1.0,Fig. 3, the traces show a long memory process, i.e., the local trend over the interval will be persistent (Korvin, 1992).

Method of successive random additions
The method of successive random additions (also known as midpoint displacement method) is a stochastic interpolation tool (for processes described by a fractional Brownian motion), which is used for generating approximately random fractals between observed data (Voss, 1988;Saupe, 1988).
The random interpolation recursive scheme allows the insertion of linearly interpolated values at the midpoint of the segment, separating the distinct points where data are given, to which a random component is added with an initial variance that decreases in every iterative or recursive level.
The initial variance is obtained from the mean square variation of the original data.This initial variance to the estimation of the mean value of the scale variations within the space gap interval between logs.The magnitude of the variance is reduced in each recursive level according to a power law determined by the Hurst coefficient (H), which is obtained for every data set.When the data values do not present a normal distribution, they are transformed into normally distributed variables before performing the interpolation process, being subsequently transformed back into their original distribution.
The random variables Z(x i ) defined at every point x i of the domain being modeled are variables that take some numerical values according to some particular probability distribution.Their spatial correlation depends (in case of translation invariance) on the vector l separating two points x i and x i + l.The set of true values z(x i ) of the variable z defining the domain being modeled is interpreted as a particular realization of the random function Z(x).
The procedure to generate a fractal distribution by applying the successive random additions method (midpoint displacement method) can be summarized by the following steps: the fractal interpolation at a given depth and between the data of two geophysical well logs is interpolated at depth h (Fig. 4).The interpolation at the points between two wells will be designed by Z i,j where i and j refer to the position and iteration order, respectively.The initial log values are Z 1,0 for well 1, and Z 2,0 for well 2. The initial variance considered in this process is given by the initial variance, σ 2 0 which is obtained from the whole data set distributed at all depths for each considered well.The process also uses the intermittence coefficient or Hurst coefficient, which is computed for each well by means of the R/S analysis technique (Hurst, 1951;Hurst et al., 1965).
The interpolation method of successive random additions is based on the fact that the incremental variance (variogram) of a random self-affine fractal trace is given by: where γ (l) is the so-called semivariogram, V H represents the variance (σ 2 ), E is the expected value of a random variable and H is the Hurst or intermittence coefficient.The step by step description of the stochastic interpolation process by the fractional Brownian motion can be summarized as follows: 1. Computation of the average initial variance (σ 2 0 ) which is characteristic of the variations between well logs.

Interpolation of the values in the midpoint interval
between wells by linear interpolation or kriging.
3. Addition of a random Gaussian number normalized at the interpolated values (or random variation) and obtained from a zero-mean normal distribution of variance, σ 2 1 where: Considering the power law scaling: for the special case of r = 1/2 we get: 4. The process is repeated recursively with all the interpolated values until the desired level of resolution is acquired.
In the n-th stage of the iteration process, the random iteration that is added to each interpolated value has the variance σ 2 n , where (Voss, 1988;Saupe, 1988): The following are the interpolation equations with R ij designating a random number drawn from a normal distribution with mean zero and unit variance.
I teration 1 I teration 4 Figures 5 and 6 show an example of fractal interpolation by the method of successive random additions for neutron porosity well logs, resulting in 17 pseudo-well logs after the fourth iteration.
A Hurst coefficient variation tracking of the shale volume trace (VSH) concerning the 17 pseudo-logs generated during the fractal interpolation process was carried out.A variation from low to high roughness between the two original wells was observed.The applied method was Rescaled Range (R/S) and the obtained results are: 0.874, 0.876, 0.876, 0.874, 0.874, 0.870, 0.866, 0.860, 0.851, 0.842, 0.834, 0.829, 0.827, 0.826, 0.827, 0.826, and 0.828.It is important to notice that low H(R/S) values represent a higher roughness of the trace.
We observed that the variation of the Hurst coefficients of the VSH traces regarding the wells distributed throughout the area of study ranged from 0.715 to 1.0 (Fig. 3), which falls within the characteristic variation range of the fGn (fractional Gaussian noise) magnitudes of the Hurst coefficient (H).For H values within the 0.5<H<1 range, a "persistent behavior"  (e.g., a positive autocorrelation) is described.For an increase occurring from time step t i−1 to t i , an increase from t i to t i+1 is very likely to happen.The Hurst exponent is also directly related to the "fractal dimension", which gives a measure of the roughness of a trace.The relationship between the fractal dimension, D, and the Hurst exponent, H, is D = 2 − H .As this equation shows, the fractal dimension is directly related to the Hurst exponent for a statistically self-similar data set.A small Hurst exponent has a higher fractal dimension and a rougher trace.A larger Hurst exponent has a smaller fractal dimension and a smoother trace.As for the variation range of the Hurst coefficient, its minimum relative values could indicate high resistivity zones and probable distributions of fluids.

Estimation and simulation of petrophysical properties
The geostatistical analysis of shale volume for several wells in the study area is presented (Tables 1 to 5).Variograms in different directions were constructed for the study of anisotropy (Gringarten andDeutsch, 1999, 2001) 1989).Several sequential Gaussian simulations for different models of shale in the cube were done.From the calculations of shale volume, geostatistical analysis can be performed to find the corresponding spatial distribution parameters.Based on these parameters, it is possible to generate estimates of its distribution in a cube using the software PETREL TM (2010).
The variogram stabilizes at the sill 0.035, which is a value close to the variance, 0.04.Another observation is that there is more continuity in the vertical direction for most of the data.The variogram model parameters are listed: nugget effect: 0.005, number of structures: 1; sill: 0.03; type of model: Gaussian, which was selected as a first smooth approximation; the ellipsoid's definition: maximum range: 520 m, midrange: 380 m, minimum range: 120 m; angles: 45, 0, 0 • (Tables 1 to 5).A mesh was constructed to estimate the fractional volume of shale in a parallelepiped.The parameters for generating the mesh are the following: number of elements in the x, y, z directions: 30, 14, 59; size of the cell in three directions: 100 m, 100 m, 10 m; minimum point coordinates: 644 500 m, 2 266 600 m, −1.925 m.
Another approach consists of generating Gaussian simulations.The parameters were as follows: number of realizations: 3; kriging type: ordinary; maximum number of conditioning values: 12; search ellipsoid definition:   1 to 5).The simulations represent different alternatives regarding the shale volume behavior.These simulations show a behavior that is more natural than the one obtained by means of ordinary kriging estimates; accordingly, they are preferred for being closer to reality.
These simulations show a correlation with the geological information and complementary geophysics that is better than the one obtained by the ordinary kriging geostatistical method.
Porosity was estimated through a sequential Gaussian simulation with a Gaussian variogram model; however, at this point, it is desirable to have a fractal approach and obtain a sequential Gaussian simulation of the effective porosity in order to reach a fractal modeling.To this end, it has to be considered that the variogram or structure function is defined by the function related to the covariance (Chilès and Delfiner, 1999).In fractal simulation, the power or fractal variogram model is used (Hewett, 1986).
Finally, in this way, an exponential variogram model was used to approach a fractal model.Then, the sequential Gaussian simulation of effective porosity, shale and permeability is carried out.As for permeability, both a porosity adjustment with measured permeability values from cores, and a sample variogram adjustment of the distribution of the permeability logarithm were used.The experimental equation, obtained from the core analysis, that was used to relate permeability and effective porosity is: Where C 1 , C 2 and C 3 are specific constants obtained for each study case.

Description of the field
The reservoir formation of Lower Paleocene to Lower Eocene age consists of turbiditic deposits grouped in three bodies (Inferior, Medium and Superior).These are limited  Based on Tyler's model (Ambrose et al., 1991), quoted in Schlumberger (2005a, b)  indicated in Fig. 7 (zone 1), where the identified discordances between some of these bodies are also shown.
Due to the fact that many of the units are complex, it is common to find different combinations of each element of this classification for a body in a given well.The definite assignation (simplified) was done according to the prevalent facies.This model allowed the performance of the spatial distribution of properties such as porosity and shaliness, which has proved very useful for the bulk analysis of each body according to the cutting values adopted for porosity and shaliness.Figures 7-11 show the porosity, shale volume and permeability distributions for each studied zone.

Results and discussion
Using the methodology detailed in the previous sections, 17 pseudo-well logs were interpolated between pairs of observed wells (Fig. 6), generated by the software PETREL TM (Petrel, 2010), obtaining three-dimensional simulations of effective porosity (PHIE), shale volume (VSH) and permeability (K) for zones 1 to 10 in the field (Figs. 7 to 11).
The properties that were used in this work for prospecting the hydrocarbon potential in the wells mentioned above were: effective porosity, shale volume and permeability.The effective porosity model  is a guide for predicting the hydrocarbon-production capacity of the reservoir.Each cell in the grid represents a value of the effective porosity in the field.The areas with yellow and red colors, which fall within 0.09-0.14,show a high level of effective porosity, while the other part of the model indicates regions with low effective porosity values.
The shale volume  represents the distribution of these petrophysical properties from the corrected version of the well log data.The grid is calibrated into fractions, which define the 3-D model in various depositional environments like the part that captures values within 20-40 % of shale content, which represents the typical regional reservoir rocks.
The permeability distribution gives another clue for the hydrocarbon potential of the field.The areas with high permeability values (yellow and red colors) within 0.1-3.0mD represent potential areas for hydrocarbon prospecting .The areas with low permeability levels allow little or no flow of hydrocarbons.
It should be noted that the results are consistent with each other, and that the regions with high effective porosity showed relatively high permeability with low shale content.In particular, zones 6, 7, 8 and 9 qualify for the presence of hydrocarbons.
On the other hand, according to the interpolation results shown in Fig. 12, it can be said that the fractal method is more powerful than the normal interpolation method, for the effective porosity distribution obtained by fractal methods models accurately the values observed in the field, whereas the distribution obtained by normal interpolation shows extreme porosity values, too high in some layers of the lateral sections and too low in the upper parts of the region under study.

Conclusions
From the results obtained in the present research, it can be seen that the sandy-clayey reservoir presents a fractal behavior as shown by the fractal analyses of the well logs.This behavior favored the application of the fractal interpolation method between neighboring wells in the field of study.This technique is used for characterizing reservoirs by means of the distribution of petrophysical properties based on well logs and core data.Through geostatistical and fractal geometry methods, predictions of the behavior of permeability, shale volume and effective porosity were obtained.
From the obtained fractal distribution, it can be said that the method of successive random additions, which was used for this purpose, is a complementary tool for the statistical characterization of reservoirs.The comparison between the fractal and normal interpolation methods justifies the fact that the fractal method is more accurate than the normal one.
The formation evaluation of the pseudo-logs obtained by fractal interpolation, using Archie's (1942) and Simandoux (1963) equations will allow the estimation of the spatial distribution of water and hydrocarbon saturations, porosity and permeability through a combination of geostatistical techniques and fractal methods.

Fig. 2 .
Fig. 2. Analysis for finding the fractal behavior of the traces used in our study, which reflected a typical fractional Gaussian noise.(a) Shale volume trace for well 2. (b) Histogram of the trace.(c) Rescaled Range analysis showing the characteristic fGn with Hurst coefficient H = 0.828.(d) Power spectrum analysis indicating power law behavior.(e) Roughness-Length method with H = 0.842.(f) Variogram analysis with Hurst coefficient H = 0.887, by using the BENOIT software.
used in our study, which reflected a typical fractional Gaussian noise.(a) Shale volume trace for well 2. (b) Histogram of the trace.(c) Rescaled Range analysis showing the characteristic fGn with Hurst coefficient H = 0.828.(d) Power spectrum analysis indicating power law behavior.(e) Roughness-Length method with H = 0.842.(f) Variogram analysis with Hurst coefficient H = 0.887, by using the BENOIT software.

Fig. 5 .Fig. 5 .
Fig. 5. Result of a fractal interpolation of neutron porosity well log data between two wells.

Fig. 6 .
Fig. 6.Map showing positions of the 14wells and the distance between them (the shortest distance between pairs of wells is 400 m.)It also indicates the 15 pseudo-wells calculated between pairs of wells.

Fig. 6 .
Fig. 6.Map showing positions of the 14 wells and the distance between them (the shortest distance between pairs of wells is 400 m.)It also indicates the 15 pseudo-wells calculated between pairs of wells.

Fig. 12 .Fig. 12 .
Fig. 12. Fractal and normal interpolation methods for the effective porosity (PHIE) distributions in the region of study.