## 1. Introduction

The microstructure of sedimentary rocks depends on their generation processes. Ninety five percent of the sedimentary rocks are clastic and carbonate formations. Clastic rocks are generated during weathering of crystalline rocks, transportation, sorting by size, etc. In contrast to the clastic formations, carbonate rocks are composed during a sedimentary deposition of organic origin or chemical precipitations. As a result of different grain shapes and posterior processes of diagenesis, carbonate formations have a complicated double-porosity microstructure. Because these sedimentary rocks are potentially hydrocarbon deposits, the determination of their porosity, saturation and permeability is required for the oil and gas field exploration and correct reserve assessment. To evaluate formation petrophysical characteristics, the interpretation of the well logging data that contain the measurements of different physical properties in boreholes is fundamental.

Geological formations, especially sedimentary rocks, can be treated as a natural composite material that consists of solid grains of different minerals and system of pores saturated with fluids (water, oil, and gas). In comparison with the synthetic-material characterization, the study of such a natural composite includes one specific additional problem that is recognition of the composite microstructure (component concentrations, shapes, spatial distributions, connectivity, size, etc) from its effective physical properties. In this chapter we present a method for characterizing the pore system structure of double porosity carbonates by using the joint inversion of well logging data. This inversion, which we have called petrophysical inversion, is based on the concept of unified microstructure model and micromechanical self-consistent method for simulating the natural-composite effective physical properties such as: elastic module, electrical conductivity, dielectric constant and thermal conductivity. We present and discuss an example of the joint nonlinear multimethod inversion of well log data for petrophysical characterization of double porosity carbonate formations.

## 2. Unified model for double – porosity medium

Generally, carbonate formations comprise two pore systems with pores of different scales, spatial distribution and shapes. According to numerous investigations (Choquette and Pray, 1970; Lucia, 1999) the carbonate porosity can be treated as a double porosity medium comprising two pore components: primary small-scale pores and secondary large-scale pores. Some examples of the carbonate rocks (Fig. 1) show existence of the complicated fracture systems generated during tectonic processes and presence of strong vug porosity. The saturating fluid flow provokes vug generation through the material dissolution in the fracture intersection or organic forms. The typical size of primary (intergranular, intercrystalline) pores does not exceed 30-50 μm while the characteristic length of secondary pores such as vugs and cracks (microfractures) are in the range of 0.5 – 50 mm (Bagrintseva, 1999). Currently the double porosity concept for modeling of rock physical properties is widely used (Barrenblatt et al., 1960; Wilson and Aifantis, 1982; Auriault and Boutin, 1994; Berryman and Wang, 2000).

The scale of the macroscopic medium description (determination of the effective macroscopic properties) corresponds to the characteristic length of well logging investigations, (0.1-0.15 m) and it is less than wavelengths of the acoustic log (0.4-0.5 m).

Since the scale of measurements on rock samples or in a borehole is much larger than the size of rock components (grains and all type pores) such a medium can be treated as a microheterogeneous material characterized by the effective physical properties. Because of the ratios between different pore scales and macroscopic scale we apply the hierarchical two step homogenization approach proposed by Auriault and Boutin (1994), and Auriault and Royer (2002).

We considered the solid frame composed by rock grains and primary pores as a matrix with effective physical parameters (non –zero conductivity and non-zero shear modulus). These matrix properties correspond to the effective properties obtained after homogenization of the heterogeneous porous frame. Then we have used this homogeneous matrix as a host in which the inclusions of different sizes and shapes representing the secondary pores are embedded. In this model we assume that all pores are completely saturated with water and are randomly distributed and oriented. Therefore the matrix and the effective double-porosity medium are treated as isotropic media. The shapes of grains, matrix and secondary pores are approximated by three-axial ellipsoids. By varying the ellipsoid aspect ratios we have modeled different secondary porosity types as cracks, vugs, and channels.

## 3. Effective medium methods

The micromechanical effective medium methods developed for composite materials have a wide application for simulation of rock properties. These methods are based on the solution

of one-particle problem for a single inclusion placed in some effective host and take into account the inclusion interactions for high concentrations. Generally two basic approaches are used for the modeling of rock geophysical parameters: the effective media approximation (EMA) and differential effective medium (DEM) methods (Budiansky and O’Connell, 1976; Berryman, 1980; Sen et al., 1981; Feng and Sen, 1985; Norris, 1985; Berge et al., 1993; Zimmerman, 1991; Hornby et al., 1994; Pribnow and Sass, 1995; Le Ravalec and Gueguen, 1996; Kazatchenko et al., 2004a, b).

The DEM method based on the sequentially repeated iterative solution of the homogenization problem works well for two-component materials in which the first constituent is accepted as a host and the second one represents inclusions placed in this host (Berryman, 1980; Norris, 1985).

However, in the case of simulating different physical properties from the unified microstructure model, the introduction of any model component as the initial host leads to serious contradictions. For example, the selection of the conductive fluid (water) saturating pore components as the host guarantees the existence of electrical conductivity for small porosities but does not allow the effective shear modulus to be calculated (Sen et al., 1981; Mendelson and Cohen, 1982; Norris et al., 1985). The host corresponding to high-resistive solid-grains inclusions (second model component) does not provide the correct simulation of the effective electrical conductivity for the low porosity range. That is why for calculating effective physical properties of geological formations we have used the symmetrical EMA method that treats equally all N components with no one distinguished as a host (Berryman, 1980, 1992; Mavko et al., 1998).

This method allows us to calculate the effective physical properties from the unified multicomponent model composed of solid grains and fluid-filled primary and secondary pores. The elements of each component are approximated by three-axis ellipsoids with different aspect ratios.

Equations for the elastic properties were obtained by Korringa et al. (1979), Berryman (1980), and Norris (1985)

where C_{i} is a volumetric concentration of the i-th component (ΣC_{i}=1), L is the elastic tensor of the effective medium, L_{i} is the elastic tensor of i-th component, and T^{(i)} is Wu’s tensor that relates the strain tensor inside an individual element of i-th component with the uniform stain field far from it (Wu, 1966).

For a medium composed of arbitrary distributed isotropic elements, the components of the tensor L_{i} are defined as_{i} and μ_{i} are the bulk and shear moduli of i-th component, δ_{ij} is Kronecker delta. In this case the system (1) can be written as:

where K and µ are the bulk and shear moduli of effective medium respectively,

For the ellipsoidal inclusion Wu´s tensor is

where I is the fourth-order isotropic identity tensor, and S^{(i)} is Eshelby’s tensor (Eshelby, 1957). The components of Eshelby’s tensor S^{(i)} for the ellipsoidal inclusions are given by the series of expressions (Eshelby, 1957). When the ellipsoid semi-axes satisfy the inequality a_{1}>a_{2}>a_{3} these expressions are:

where υ is Poisson’s coefficient of effective medium, J_{j}, J_{ik} are integrals given by

In these expressions F(θ,k) and E(θ,k) are the elliptic integrals of the first and second kind where

We used Carlson's algorithm (Carlson, 1979) to calculate the elliptical integrals F(θ,k) and E(θ,k). The components of the Wu’s tensor T(i) (Wu, 1966) were calculated numerically using equations (4).

Electrical conductivity equations for isotropic medium with ellipsoidal components were obtained by Sen et al. (1981) and Norris et al. (1985)

where_{i} and σ are electrical conductivities of the i-th component and effective medium, respectively. The tensor D^{(i)} is analog of the Wu’s tensor for the electromagnetic field.

where n^{(i)} is depolarization factor of the i-component (Stratton, 1941). The depolarization factor has only three none zero components n_{kk}^{(i)}. Using the appropriate ellipsoidal coordinate system (Landau and Lifshitz, 1960) these components can be determined as following (index i of the constituents is omitted)

and calculated using the elliptic integrals F(θ,k) and E(θ,k)

here a_{1}, a_{2,} and a_{3} are the ellipsoid semi-axis.

The thermal conductivity and dielectrical permittivity for composite medium are described by equation (5) where σ_{i} and σ are substituted by λ_{i}, ε_{i} and by λ, ε, respectively where λ_{i}, ε_{i} are thermal conductivity and dielectrical permittivity of i-th component and λ, ε are thermal conductivity and dielectrical permittivity of effective medium

Equations (2), (5), and (9) are solved numerically by iterations.

## 4. Simulation of physical properties

By using the unified model and EM method we have modeled the following effective properties of a double porosity medium: acoustic P- and S-wave velocities, electrical conductivity, dielectric permittivity, and thermal conductivity as functions of the primary and secondary porosity, and shapes of secondary pores (aspect ratios). In this section we consider a matrix model and determination of its parameters, and then modeling of the physical properties of a double-porosity medium.

### 4.1. Matrix model

As it was mentioned, for simulating physical properties of the double porosity formation we apply the hierarchical two step homogenization. At the first step the matrix effective properties are obtained. The matrix is presented by a two-component composite medium of solid grains and pores. The grains of high resistivity form the solid frame. The pores are saturated by a conductive fluid. The elements of each component are approximated by three-axis ellipsoids with different aspect ratios which are introduced as functions of porosity.

The introduction of the varied aspect ratios is based on the supposition that the microstructures of porous rocks (pore-size distribution and its connection network) for high and low porosities are different. The variations of the aspect ratios can be obtained from two classes of experimental data. The first includes data about the porosity and microstructural statistical characteristics measured on cores. The second class (used in this work) presents the acoustic velocities and electrical conductivity from cores or well log data. In this case the aspect ratios are determined by adjusting the calculated parameters (acoustic velocities and electrical conductivity) with experimental data. For acoustic-wave velocities (V_{P} and V_{S}) we applied the generalized regression equation for a limestone formation saturated with water for porosities from 0.01 to 0.32 (Mavko et al., 1998).

where V_{Pm} and V_{Sm} are matrix P- and S-wave velocities (km/s), ρ_{m} is matrix bulk density (g/cm^{3}),_{f} is the bulk density of a saturating fluid (that in the considered case is equal to 1 g/cm^{3}).

The relationship of the electrical conductivity to a primary porosity for the carbonate formations that has no secondary porosity as microfractures and vugs, was obtained by Kazatchenko and Mousatov (2002) based on the statistical analysis of numerous experimental data. The formation factor for this type of rock is described by Archie’s law (Archie, 1942)

with the cementation exponent m=2, here σ^{*} and σ_{f} are electrical conductivities of a medium and a saturating fluid respectively.

Using equations (2) - (8) for the components of corresponding tensors obtained for ellipsoids, we calculated the P-wave velocity and conductivity as functions of porosity. The minimization of a fitting error between experimental and synthetic data allows determination of the pore and grain geometry (i.e., values of the aspect ratios) over the porosity range of 2-31%. These aspect ratios were used for calculating S-wave velocities and dielectric permittivity to compare them with the regression equation (10) and (12) to further confirm the validity of the model structure obtained. Semi-empirical equation (12) (CRIM) was obtained for homogeneous rocks with one-pore system (Guéguen and Palciauskas 1994)

Saturating water parameters are: ρ_{f}=1 g/cm^{3}, K_{f}=2.25 GPa, electrical conductivity σ_{f}=1 Ωm, and the ratio of the water and grain conductivities is taken to be 105. The limestone and water dielectrical permittivity are ε_{l} =7.5 and ε_{w}=80+i0.05.

First, we assumed constant aspect ratios of pore and grain elements on the all porosity interval (2-31%). For this situation, the best coincidence of theoretical and experimental data is achieved when the aspect ratios of grains and pores are α_{1g}=0.3048, α_{2g}=0.0903 and α_{1p}=0.0803, α_{2p}=0.0252. The matching of the acoustic data is satisfactory, with a relative error less than 2%. However, the conductivity relative error is rather high and corresponds to 30% in the porosity interval 5-10%. For porosities lower than 5% the predicted conductivity is reduced considerably and tends to the percolation threshold.

The variation of the aspect ratios allows the predicted and experimental data in the porosity interval 3-31% to be adjusted with a high degree of accuracy (Fig. 2 A and B). The relative error for the P-wave velocity is less than 0.05%, and for the conductivity does not exceed 5% for porosities over the range of 4-31%.

The EMA using spheroids with non variable aspect ratios was rarely applied to model electrical properties, because it predicted vanishing conductivity (percolation threshold) at relatively high porosity values. Allowing variation of the pore and grain aspect ratios moves the percolation threshold to lower porosities near 2% in agreement with measured data for carbonate collections.

The calculated S-wave velocities coincide with experimental equation both for constant and for variable aspect ratios (Fig. 2C). The maximum relative errors are 2% and 3% respectively. The dielectrical permittivities obtained with variable shapes of pores and grains are in a good agreement with CRIM equation. Small differences between them increases with a porosity value however it does not exceed 7% for porosity 30% (Fig. 2 D and E).

For high porosities (28-31%) the pores have ellipsoidal form with aspect ratios α_{1p}=0.300 and α_{2p}=0.097. When the porosity decreases from 31% to 2% the ellipsoids transform into the needle-form pores with the small aspect ratios α_{1p}=0.065 and α_{2p}=0.022. This pore geometry maintains the electrical conductivity for low concentrations of conductive component. The solid grain geometry changes from the ellipsoids with α_{1g}=0.408, α_{2g}=0.269, to penny-shapes characterized by α_{1g}=0.984, α_{2g}=0.098. The variations of aspect ratios with porosity are approximated by polynomials of the second power with the correlation coefficient

The predicted variations of shapes for skeleton grains and saturated pores correspond to the recently reported experimental data for carbonate formations on pore microstructure changes depending on porosity (Song et al., 2002).

For the matrix thermal conductivity we have applied the regression equation from Lichtenecker’s model (Pribnow and Sass, 1995; Schoen, 1996)

where λ_{m} and λ_{f} are thermal conductivities of matrix and saturating water, λ_{f} =0.63 W/m*K.

### 4.2. Model of double-porosity rock

On the second step of the hierarchical homogenization we considered the solid frame composed by rock grains and primary pores as a matrix with effective physical parameters (non –zero conductivity and non-zero shear modulus). The secondary pores correspond to inclusions placed into the matrix and they are approximated by three-axial ellipsoids with aspect ratios α_{1} and α_{2}. Selection of the aspect ratios allows us to model different type of secondary porosity such as vugs (close to sphere inclusions), vugs connected by channels (quasi needle shapes), and cracks (flattened ellipsoids). The shape of the matrix component (taking into account that the matrix is an isotropic homogeneous medium) is given by spheres.

We have simulated the compressional and shear wave velocities, electrical and thermal conductivities for two cases: (1) the secondary porosity comprised by one type of pore-shapes (unimodal secondary-porosity shapes) and (2) the secondary porosity consisted of two components with different aspect ratios (bimodal secondary-porosity shapes). The effective physical properties were calculated for different inclusion aspect ratios α_{1,2}=1-0.001, and values of the matrix and secondary porosities: ϕ_{m}=0.02-0.24 and ϕ_{s}=0.0025-0.08. The dielectrical permittivity was simulated only for the first case.

Because the experimental measurements in double porosity rocks provide the value of total porosity, for quantitative evaluation of the secondary-pore influence it is convenient to calculate deviations between the effective physical properties for two media with the same total porosity ϕ_{t1}=ϕ_{t2} but different concentrations of the primary and secondary pores. The first medium has only matrix porosity ϕ_{m} (ϕ_{t1}=ϕ_{m1}), while another one has both the matrix (ϕ_{m2}) and secondary (ϕ_{s}) porosities. The total porosity of the second medium is ϕ_{t2}=ϕ_{m2}+ϕ_{s}- ϕ_{m2}ϕ_{s}. The deviations of the acoustic velocities ΔV_{P} and ΔV_{S}, electrical Δσ, thermal conductivities Δλ, and dielectrical permittivity Δε are ΔV_{P}=V_{P}-V_{Pm}, ΔV_{S}=V_{S}-V_{Sm}, Δσ=ln(σ)-ln(σ_{m}), Δλ=λ-λ_{m}, Δε=ε-ε_{m}, where V_{P}, V_{S}, σ, λ, and ε are the effective parameters of double porosity medium with the total porosity ϕ_{t2}, V_{Pm}, V_{Sm}, σ_{m}, λ_{m}, and ε_{m} are parameters of a medium with the total porosity ϕ_{t1}.

It should be noted that instead of the fracture density (

### 4.3. Secondary-porosity with one type of pore-shapes

We have calculated the differences of the effective elastic, electrical, thermal and dielectrical parameters for secondary porosity ϕ_{s}=0.02 and two matrix porosities ϕ_{m}=0.08, 0.18 (ϕ_{t2}=0.1, 0.2). Depending on the aspect ratios, the parameters studied lie below (negative deviations) or above (positive deviations) the corresponding matrix regression equations (Fig. 3). The electrical conductivity has considerable sensitivity to the secondary pore shape and value. For the matrix porosity ϕ_{m}=0.08 the absolute deviation Δσ changes from +0.5 for cracks (α_{1}=1, α_{2}=0.001) to -0.3 for vugs (α_{1}= α_{2}=1) that corresponds to the relative sensitivity of +12% and -8% (Fig. 3A). Both cracks and vugs affect the effective electrical conductivity approximately in the same way. However, the relative sensitivity of the electrical conductivity drops up to ±5% when the matrix porosity increases.

Cracks increase and vugs decrease the dielectric permittivity (Fig. 3 B). Its influence depends directly on the value of the secondary pores. Δε may achieve 5-6% for ϕ_{s}=0.02 and matrix porosity ϕ_{m}=0.08. However, for ϕ_{m}=0.18 the difference Δε is less than 2%.

The thermal conductivity demonstrates a very low relation with the secondary-pore shapes (Fig. 3 C). The difference Δλ does not exceed 1% even for model with the ϕ_{s}=0.02 and ϕ_{m}=0.08 and the parameter Δλ decreases with increasing matrix porosity.

Presence of the secondary pores with α_{1}<α_{2} results in a decrease of the acoustic velocities with respect to the matrix velocities (Fig. 3 D and E). Strongly flattened ellipsoids with α_{1}<<α_{2} (cracks) cause considerable relative deviations of velocities: 30% and 85% for P- wave and S-wave respectively. The sensitivity of the S-waves to cracks is about three times bigger than that of the P waves. The secondary inclusions with α_{1}≈α_{2} slightly increase the effective velocities (in the range of 1%) in comparison with the matrix values. If the secondary porosity remains the same, the increment of the matrix porosity leads to smoothing the secondary porosity effect.

Variations of the ΔV_{P}, Δσ, and their zero-isolines versus the aspect ratios of the secondary pore are plotted in Fig. 4. We confined our examination of velocity variation to ΔV_{P} parameter because ΔV_{S} versus aspect ratios has the same character as ΔV_{P}. The Δλ and Δε have a small sensitivities to the secondary-pore shapes. When the secondary ellipsoids have aspect ratios α_{1} and α_{2} bigger than 0.03 the effective electrical conductivity is lower than the matrix conductivity (Δσ<0). Outside this region (Fig. 4), when aspect ratios decrease, i.e. the secondary pores become elongated in one or two dimensions, the electrical conductivity exceeds the matrix one (Δσ>0).

The central area on the plots for P-wave velocity (Fig. 4), which correspond to nearly sphere or prolate spheroids, is characterized by velocities bigger than the matrix one. A medium with such inclusion shapes is more rigid than the matrix with the porosity ϕ_{t1}= ϕ_{m1} (ΔV_{P}>0). In two periphery zones, where the secondary-pore shapes approach flattened ellipsoids, there is a reduction of the effective compressional wave velocities with respect to the matrix.

Presence of the secondary pores corresponded to the central zone in Fig 4 generally increases the thermal conductivity. However, strongly flattened ellipsoids (α_{1}=0.03-1, α_{2}=0.001-0.03 and α_{1}=0.001-0.03, α_{2}=0.03-1) make the thermal conductivity lower than that of the matrix.

It should be noted that for certain combinations of the α_{1}, α_{2}, ϕ_{m}, and ϕ_{s} the effective parameters of double porosity medium is equal to the matrix parameters (the zero-isolines in Fig. 4). The distribution of these isolines has a quasi orthogonal behavior and, consequently, the responses of the acoustic and electrical parameters can be considered independent from each one. The elongation of inclusions determines the electrical conductivity changes. Elongated inclusions improve the electrical current propagation but do not affect the elastic module. On the other hand, the elastic velocities are sensitive to the flatness of inclusions (two dimensional elongations). This secondary-porosity type corresponds to cracks or microfractures that make the medium softer and result in significant decrease of the P- and S-wave velocities.

### 4.4. Classification of the secondary-porosity types

Based on the ΔV_{P} and Δσ responses we have divided the inclusion shapes into four general groups that are associated with the specific types of secondary porosity. Due to the weak inclusion effect on the thermal conductivity and dielectric permittivity we do not consider these parameters in the secondary-pore classification. We used the zero-isolines (boundaries where the signs of the ΔV_{p} and Δσ deviations are changed) for delimiting the pore-shape groups (Fig 4). Such a pore-shape separation is convenient because the position of sign-change boundaries weakly depends on the values of matrix and secondary porosities.

The first group (I) of inclusions is announced to Δσ<0 and ΔV_{P}>0. It corresponds to the vuggy porosity. Inclusions from the second group (II), where Δσ>0 and ΔV_{P}>0 have needle shapes and can be considered as channel porosity (Fig. 4). Flattened ellipsoids with aspect ratios in the area III (Δσ>0 and ΔV_{P}<0) describe the crack porosity type. Fourth group (IV) is characterized by Δσ<0 and ΔV_{P}<0. This pore–shape group forms a transition porosity type between vuggy and crack porosities. We call this quasi-vuggy porosity type.

The deviations Δσ, Δλ, ΔV_{P}, and Δε as functions of the total and secondary porosity for each porosity types are shown in Fig. 5-9. These graphs were calculated for the following aspect ratios corresponding to vugs (α_{1}=α_{2}=1), channels (α_{1}=α_{2}=0.01), cracks (α_{1}=1, α_{2}=0.001), and quasi-vugs (α_{1}=1, α_{2}=0.1). For all pore-shapes considered, the deviations of the physical properties are proportional to the secondary porosity and inversely depend on the matrix porosity. In addition to the general tendency, each parameter has specific features:

Electrical conductivity (Fig. 5). As noted earlier, electrical conductivity has a high sensitivity to the presence of secondary pores for all porosity types. However, the sensitivity strongly decreases with increasing matrix porosity. For matrix-porosity value above 0.1-0.15, the effective conductivity rapidly approaches the matrix conductivity and practically coincides with it for ϕ

_{m}>0.2.Thermal conductivity (Fig. 6). For all porosity types the deviations Δλ are low and may be neglected in the range of matrix and secondary porosity values considered here.

Dielectric permittivity has complicated relations with total and secondary porosities for models with channels and quasi-vugs. These relations result in the increase of the dielectric permittivity for high total porosity (effect of cracks) and the decrease of the dielectric permittivity for low total porosity (effect of vugs) (Fig. 7 B and D).

Acoustic wave velocities (Fig. 8 and 9). The acoustic wave velocities have variable sensitivities to the types of secondary porosity. In the cases of vug and channel porosities (Fig. 8 A, B and Fig. 9 A, B) the acoustic velocities are weakly depended on the matrix porosity.

It should be noted that a direct relation between the ΔV_{P} and the matrix porosity occur for channel porosity. This fact can be explained by decreasing the effective V_{P} more slowly than the matrix V_{Pm} with the growth of the matrix porosity. Cracks and quasi-vugs strongly affect the acoustic wave velocities. Small crack concentration (ϕ_{s}=0.005) leads to significant deviations of the effective velocities from the matrix. The sensitivities of velocities to cracks and quasi-vugs are maintained high in the whole matrix-porosity interval despite of weakly decreasing with matrix porosity increase.

### 4.5. Secondary porosity composed by the mixture of two pore-shape types

To calculate the effective properties of media with two types of secondary pores we applied the EMA scheme for three-component: matrix represented by elements of spherical shape and secondary pores formed two different groups of ellipsoids.

Small concentrations of vugs in a medium with channel porosity (Fig. 10) or channels in a rock with vuggy porosity (Fig. 11) keep Δσ, Δλ, and ΔV_{P} close to the deviations corresponding to the dominant porosity type. Presence of vugs reduces the electrical conductivity of rocks with the channel-porosity type and, conversely, introducing channels in vuggy-porosity media enhances electrical conductivity. When the values of channel and vuggy component are the same, the electrical conductivity coincides with the matrix regression equation (

Addition of a low concentration of cracks in media with dominant vuggy or channel porosity leads to significant changes of the acoustic velocities (Fig. 12 B and 13 B). For example, a crack fraction of 0.005 added to a rock with the 0.08 porosity of vugs or channels is enough to reduce the effective velocities to values lower than the matrix velocities. In this situation the velocity response corresponds to the cracked medium and does not reflect the vuggy or channel porosities. Presence of crack weakly increases the electrical conductivity of the channel or vuggy porosity media (Fig. 12 A and 13 A) and it does not affect the general behavior of Δσ-curves defined by their dominant porosity types. Furthermore, the influence of crack on the conductivity decreases for increasing vuggy or channel porosities. The thermal conductivity response depends mostly on the value of secondary porosity (vuggy or channel) and does not vary when cracks are added.

Thus a secondary porosity with two pore-shape components that may occur in real carbonate formations produces specific variations in the effective physical properties and significantly changes the relationship between the electrical and acoustic parameters obtained. When vug porosity is equal to the channel porosity, effective electrical conductivity coincides with the matrix conductivity. In this case, only acoustic velocities indicate the presence of secondary pores (ΔV_{P}>0). Addition of a low crack concentration in vuggy- or channel-porosity rocks changes completely the velocity response, resulting in characterization of such media as fractured only. In this case electrical conductivity response will indicate the dominant type of secondary pores. A specific combination of secondary-porosity types can lead to the formal equality of one effective parameter to the matrix value, while the response of the other physical property will show the dominant porosity type. Therefore the joint analysis of the acoustic wave velocities and electrical conductivity allow us to estimate correctly the secondary- porosity with one or two types of pore-shapes.

The calculation results for the models completely saturated by water, show different sensitivities of the acoustic-wave velocities and electrical conductivity to the inclusion shapes, matrix and secondary porosities. The electrical conductivity demonstrates the similar sensitivity to both the vuggy and crack porosity types. While the acoustic wave velocities are significantly affected by flattened inclusions and have very low sensitivity to nearly spherical inclusions. The thermal conductivity is characterized by a low sensitivity to the water-saturated secondary pores for all considered shapes.

Based on the results obtained we propose a new secondary-porosity classification: vugs (close to sphere inclusions), quasi-vugs (weakly oblate inclusions), channels (needle-shape inclusions), and cracks (strongly flattened inclusions).

The variations of effective characteristics when the secondary porosity was constituted by mixtures of two pore-shape types (such as vugs and cracks, channels and cracks, vugs and channels) are also assessed. The joint analysis of the velocities and conductivities allows us to determine correctly the composition of the secondary porosity constituted by two pore-shape components.

## 5. Joint inversion for the rock microstructure determination

In this work, based on the unified microstructure model, we perform the joint nonlinear multimethod inversion of well log data. A set of experimental data includes measurements of different physical properties of geological formations (total porosity, density, P- and S-wave velocities, electrical conductivity, and natural radioactivity) carried out along a borehole with high density of observations.

The inversion procedure consists in minimizing a cost function that includes the sum of weighted square differences between the measured well-logging data and theoretical logs calculated using unified model and EMA method. Taking into account that this optimization procedure is an ill-posed problem because the experimental data contain noise related to equipment accuracy and environmental conditions, we have included additional regularization functional containing complementary information (initial model parameters and ranges of their variations) in the cost function. The petrophysical parameters determined by the inversion are the mineral concentrations, matrix porosity, secondary porosity and secondary porosity type (aspect ratios of spheroids approximating secondary pores). The physical properties of mineral components and fluid are assumed as known and can be adjusted by the posterior analysis of fitting-error distributions for each log.

To determine the pore microstructure parameters we have formulated the joint inversion procedure of well log data as solving an optimization problem where the following quadratic cost function is minimized for each point of measurements

here the vectors d_{obs} and d(m) are the measured and theoretically calculated well log data, respectively. These vectors can be presented in the form

where their components correspond to the logarithmically scaled conventional well logs used in the inversion: P-wave transit time (the parameter is inverse to velocity) – t_{P}, S-wave transit time – t_{S}, resistivity measured by MSFL tool - r_{MSFL}, resistivity measured by LLD tool - r_{LLD}, neutron porosity - ϕ_{t}, density - ρ, gamma ray - γ, and the log of photoelectric effect – PEF. We have found that the use of the logarithmic scale is more suitable for the joint inversion of different logs because they are positive and some of them are characterized by

large variations (resistivity, porosity, gamma values) and close to lognormal distributions (or at least non-Gaussian distributions). In this case the datum residual term of equation (16) represents normalized deviations between calculated and measured data and allows us to avoid the problem related to scales and units of different measured quantities.

The variable m is the vector of unknown pore microstructure parameters which have to be found by minimizing the cost function F(z_{i})

where ϕ_{m} is the matrix porosity, ϕ_{s1} is the secondary porosity of the first type (microfractures), ϕ_{s2} is the secondary porosity of the second type (vugs and channels), α_{s1} and α_{s2} are the aspect ratios of spheroids (pore shape) of the first and second types of the secondary porosity respectively, V_{sh} and V_{d} are the shale and dolomite volumes respectively. The vector m_{0} represents a reference model that can be given using a priori information

W_{d} is a diagonal matrix of weight coefficients which allow us to account different error scales and distributions for each log. These coefficients are a nonlinear functions of the relationships between the physical characteristics (given logs) and model parameters m, and accuracy of measurements for each log. A diagonal matrix W_{m} assigns the weights for model parameters based on a priori information. We adjust the component values of both matrices using analysis of the posterior error distributions for each log.

The scalar factor λ is a regularization parameter that introduces a relative weight between the misfit term (first term) and Tikhonov’s stabilizer (second term) of the cost function (Tikhonov and Arsenin, 1977). The optimal value of this parameter has to provide the minimal deviation from the reference model m_{0} (maximal influence of a priori information) and keeps the misfit between the calculated and measured logs within a prescribed error values (Zhdanov, 2002).

To solve the optimization problem for the cost function (equation 16) we have applied the Levenberg-Marquardt method (Levenberg, 1944) that provides the stable global minimum search without calculating functional derivatives. This method is especially effective when the minimized functional has a complex surface as in the case of joint inversion of different physical properties.

We have tested the presented joint inversion technique on the theoretically simulated well logs with prescribed levels of noise for synthetic models (matrix and secondary porosities, secondary pore shapes) and on experimental data by comparing inversion results with geological and petrophysical information (core data, NMR, FMI).

### 5.1. Synthetic model

Here we present a verification of the inversion convergence on a theoretical model that consists of grains, matrix pores, cracks and vugs approximated by spheroidal inclusions. The parameters of the solid grains composing matrix (density, P- and S-wave velocities) correspond to limestone. The matrix and secondary pores are saturated with water. The

matrix, crack, and vuggy porosities were changed in the ranges of 3 - 5%, 0.1 - 1.5%, and 3-6%, respectively. The crack and vug aspect ratios were varied from 10^{-3} up to 5*10^{-3}, and 0.1-0.7, respectively. The shapes of matrix pores and grains are functions of the matrix porosity (equation 13). The values of all parameters were chosen in the mentioned intervals using uniform random distribution (the Monte-Carlo method). This model have been generated by 1000 realizations. For each realization the P- and S-wave transit times, electrical resistivity, density, and total porosity were calculated adding a random and normally distributed noise of 3% for probability 0.97 which correspond to well logging measuring errors. Then we have applied the joint inversion of the simulated logs to reconstruct the model parameters. The comparison of the matrix, vuggy and fracture porosity, and aspect ratios of vugs and fractures obtained is presented in figures 14. The dotted lines mark the prescribe error-intervals: absolute error ±0.003 for porosities and relative error ±10% for aspect ratios. More than 50% of the matrix and secondary porosity determined by inversion are in the range of ±0.002 and the fracture and vug aspect ratios are in the range of ±10%.

### 5.2. Model of carbonate formation with thin shale layers

The volume of shale in the carbonate formations is not high generally and it does not make difficult the potential zone detection using logging data. On the other hand the presence of shale can lead to erroneous results in the determination of secondary pore types if the model of pure carbonate is applied for the joint inversion of logs. In carbonate sediments the shale component are mostly presented as thin layers and in some special facies, it can be arbitrary distributed as inclusions in the matrix or in the secondary pores (Scholle and Ulmer-Scholle, 2003).

Based on this geological data we have introduced a model for carbonate formations containing shale that corresponds to a transversely isotropic medium composed of intercalated layers of pure carbonates and shale. In this model the axis of macroscopic anisotropy (appeared due to formation texture) is perpendicular to the stratification plan and the effective physical properties are described by tensors. For the vertical borehole and horizontal layers the tensors of effective parameters have only two principal components which we will name as vertical (perpendicular to layers) and horizontal (parallel to layers) components.

Taking into account that in the acoustic logs the vertical velocities of elastic wave are measured, we can apply the average time equation for the P- and S – waves’ travel times

where t_{υ}^{P}, t_{c}^{P}, t_{sh}^{P}, t_{υ}^{S}, t_{c}^{S}, and t_{sh}^{S} are the travel times of vertical effective velocity component, travel times in carbonates and shales for P - and S – waves, respectively, and V_{sh} is the shale volume. The resistivity well logs in thick layers due to anisotropy paradox measure the horizontal component of conductivity tensor

where r_{h}, r_{c} and r_{sh} are the horizontal component of the effective resistivity tensor, carbonate and shale resistivities, respectively.

The effective density ρ_{υ} and total porosity ϕ_{υ} in this model can be given as follows

where ρ_{c}, ρ_{sh}, ϕ_{c} and ϕ_{sh} are the densities and total porosities of carbonates and shales, respectively. The relationships of the gamma ray and photoelectric (PEF) logs as functions of the shale volume was approximated by a linear commonly used equations (Bassiouni, 1994)

where γ_{υ}, γ_{sh}, and γ_{c}, are the effective gamma, and gamma values for the shale and clean carbonates, respectively; PEF_{υ}, PEF_{sh}, and PEF_{c} are the effective PEF parameter, and PEF parametres for the shale and clean carbonates, respectively.

The effective elastic moduli (P – and S – waves’ velocities) and electrical resistivity for the pure carbonate component are found by using the modeling approach described above for the simulation of physical properties starting from the unified pore-space model of double porosity medium. The shale physical properties as well mineral and fluid characteristics can be considered as known and then they are adjusted by the posterior analysis of inversion error distributions for each log.

## 6. Examples of the well log inversion

Below we present the results of joint inversion for two boreholes in carbonate formations characterized by different types of secondary porosity.

The first formation is described as compact limestone that are characterized by the low neutron porosity and composed of matrix porosity and microfractures with porosity in the range of 0.01-0.03 (Fig. 15). The matrix porosity does not exceed 0.01-0.02. The average fitting errors between each calculated and measured log is about 5%. The interval x150-x160 m contains microfractures with the porosity about 0.015-0.02 and the relativity high aspect ratios 0.01-0.05 that can be associated with open microfractures affected by dissolution. In this interval the presence of gas was detected. The zone x160-x180 m corresponds to the approved oil production interval with the secondary porosity up to 0.04-0.06 formed by open microfractures and vugs. The vug shapes represents oblate spheroids with the aspect ratios in the range from 0.01 up to 0.1. The results of inversion indicate that the production interval is associated with secondary-porosity increasing and can be extended downwards up to depth of x185 m. The additional perspective zone can be marked in the interval of x245-x265 m. In the core sampling interval x210-x218 m the pore microstructure characteristics obtained by the inversion correspond to compact formations (the matrix porosity is lower than 0.01) with the secondary porosity of 0.01 represented mostly by cracks. This result coincides with the core description as compact fractured mudstone. The second inversion example represents high porosity carbonate formations with the vuggy porosity type (Fig. 16). This formation is composed of dolomite and limestone with low enough shale volume. At the depth of core sampling (x370 m) the dolomite concentration determined by inversion coincides with one obtained from core. The matrix porosity varies in the range of 0.02-0.04 while the secondary vuggy porosity changes from 0.02 up to 0.16.

The misfit distributions demonstrate that the average errors between the measured and calculated logs (except the gamma log) are about 5%. High dispersion of the corrected gamma log is related to the incomplete correction of uranium concentration and probable presence of disperse clay that is not took into account in the applied model. The inversion results show a good agreement with the images obtained by FMI tool (Fig. 16). The results of joint inversion obtained for various boreholes from vuggy and fractured carbonate reservoirs show a good correspondence with core data, image log, and geological descriptions. The output inversion parameters as the matrix porosity, secondary porosity, and aspect ratios can be considered as important quantitative characteristics of pore microstructure to improve classification of carbonate lithotypes, processes of secondary-pore generation (fracturing, dolomitization, dissolution), and permeability prediction.

## 7. Conclusions

We have presented a technique for joint simulating the acoustic velocities, electrical and thermal conductivities, and dielectric permittivity of double porosity rocks treated as a natural composite material. The technique is based on the unified microstructure model and hierarchical two step homogenization using the EMA symmetrical approach. The peculiarity of the model consists in calculating the effective properties of the porous matrix with the variable aspect ratios for ellipsoids approximated by grains and pores which depend on the porosity value (first step of homogenization). This approach provides simultaneously a non-zero shear modulus for high porosity and electrical conductivity for low porosity and allows us to simulate the effective matrix properties with a high accuracy.

In the second homogenization step we calculate the effective properties for the double porosity medium where the secondary pores of different shapes are embedded into the effective matrix. The presence of secondary pores changes significantly the acoustic and electrical properties of the double porosity medium in comparison with composites containing just one pore system despite the same total porosity values. For correct estimation of the secondary-porosity with one or two types of pore-shapes, both acoustic wave velocities and electrical conductivity have to be used.

Based on the modeling results we have developed the method of nonlinear joint inversion of measured physical properties for determining the microstructure parameters of natural composite porous material. The method has been successfully applied for the integral interpretation of experimental well log data for petrophysical characterization of double porosity carbonate formations containing cracks and vugs.

The matrix and secondary porosities, and shapes of secondary pores obtained as the results of joint inversion of well logs can be considered as important quantitative characteristics of the pore microstructure to improve evaluation of carbonate formations and allow us (1) to characterize carbonate lithotypes and processes of secondary-pore generation (fracturing, dolomitization, dissolution), (2) to determine saturation of double-porosity rocks, (3) to determine geomechanical parameters and calibrate seismic data, (4) to assess probability of the secondary pore percolation for predicting formation permeability.