Creating A Digital Pore Scale Rock Modeling Biology Essay

Published: Last Edited:

This essay has been submitted by a student. This is not an example of the work written by our professional essay writers.

In this section, the theoretical findings of the authors Øren et al. [12, 13, 14, 15], Pilotti [23], Voronoi [13] and Hilfer [25] are reviewed and presented.

The necessary input data for creating a digital pore scale rock modeling is mainly obtained from thin sections. After applying image analysis techniques, the actual description of the pore space of sandstone rocks is obtained. However, for the purpose of this work the input data will be generated from a statistical function (random function) of well known outcrops type of rocks in order to create stochastic models.

The basic idea is to emulate the real sandstone rock and its main forming geological processes, such as sedimentation, compaction and diagenesis [12].

[Two lines were erased and one rephrased]

2.1.1. Sandstone geological processes

Three processes describe the main forming steps for the sediments to form a sandstone rock.


From the grain size distribution curve, grains are chosen randomly depending on their radius and at the same time a numbering process starts. The grains selected are dropped into a box making the grain bed.

The procedure for modeling ellipsoidal grains is as follows;

From the grain size distribution curve, a grain radius is picked.

The grain sphericity is assigned randomly. (is a user defined value)

Once a grain is about to fall onto the box, all the radii of the existing grains increase the radii by an mount limited by the aspect ratio assigned in the previous step.

When the grain is dropped, this one will reduce to a point and will fall onto the grainbed until its find its stable or minimum global position.

Reset all the grains radii to its initial value and start all over again for the next one.

It is assumed that all the ellipsoidal grains could be subjected to either low energy environment or high energy environment. Also, they are organized to be with their major section horizontal.


Compaction is obtained from the overburden (vertical direction) and from tectonics movements coming from lateral forces. Once this process starts, the volume of the matrix will eventually decreased and would be considered as porosity loss forcing the expulsion of the fluids in the pores [65].

From the center of the grains, the vertical coordinate is changed towards the centre x-y plane of the model [56]. The rearrangement and compaction effect can be seen in the following equation:

where z is the desirable upright position, , is the compaction factor, Zo is the original vertical co-ordinate, and is a variable that emulates how the grain rearranged. More information about the compaction modeling can be found elsewhere [12, 14].


Diagenesis is considered to start at the moment that sediment is deposited and buried by other sediments covering all the physical, chemical and biological processes that modify the sediments during the burial process. It continues until high pressure and temperature drastically change the structure and mineralogical composition of the rock, which is referred to as metamorphism.

In this work, some of the diagenetic processes are modeled, such as dissolution of feldspars, quartz cement overgrowth, carbonate, mica and clay growth. There are two methods used to model cement overgrowth, Pillotti's technique [23] and Voronoi [13]. The Pilloti technique was used in most of the geological model whereas Voronoi in just a few.

This technique inserts the original grain (sphere) inside an ellipsoid and performs a user-defined tangent planes cutting procedure to the original sphere. The tangent planes govern the angularity of the grains (see figure 2.1). The goal is to simulate the cement overgrowth that commonly occurs during diagenesis process of sandstone rocks.

The main goal is achieved applying the following procedure:

Target porosity or the desire porosity of the model has to be set and defined. This number will be the most dominant parameter during the iteration procedure.

The user has to define the way the tangent planes are implemented in the model. The options are fixed and randomly distributed in a range. In addition, the inner radius of the ellipsoid has to be determined through the parameter sphericity factor. This parameter is as well user-defined and is given by the ratio between the inner radius of the ellipsoid and the radius of the original grain.

Once the previous steps are completed, the process starts simulating the cement overgrowth through the surface expansion of each grain that belongs to the porous medium. The end of the algorithm is reached when the value of porosity from the modified porous medium is computed and compared with the desire porosity value defined in step i. If is necessary, the procedure is iterated again.


Figure 2.1. a) Representation of a 2D schematic of the modeling of non-spherical grains using Pilloti technique. b) Translation from cutting lines to cutting planes.

In contrast, using Schwartz [24] algorithm for quartz cement overgrowth, Voronoi's technique starts the process.

Voronoi cells are polyhedron carrying on the volume closest to each grain centre in the porous medium. The equation 2.2 shows how the new Voronoi radius is found:

where Ro is the deposited grain radius, l(r) is the distance between the surface originated by the Voronoi polyhedron along the direction of r and the original grain surface, is a constant that controls the amount of cement growth and known as the growth exponent that controls the direction of growth.

Each of the previously defined parameters in equation2.2 plays an important role in Voronoi's technique but the growth exponent could certainly give a clearly insight of the model. When have positive values, the growth of the quartz cement is shown in the direction of the pore bodies or large l(r) values whereas negative values will influence the growth through the pore throats or small l(r) values. When the pore throat space is not affected (positive ), the interconnection of the pores are available even al low porosities. In contrast, when the pore bodies space are not affected (negative), the pore throats tends to be clogged as the porosity decrease and therefore the tortuosity of the porous medium will increase. Finally, if is zero, then the growth occurs equally in all directions i.e. both pore bodies and pore throats.

Other important effect during diagenesis is the diagenetic clays. The morphology of the common diagenetic clays can be divided into three broad categories: pore-filling, pore-lining and pore-bridging.

The pore-filling and pore-lining are modeled in a similar way. Both are controlled by the clay clustering factor. When this factor is high, the pore-filling clays are more likely to be modeled in the form of hexagonal booklets of kaolinite that precipitates in randomly pore bodies that already contains clays. Whereas when the factor is low, the pore-lining clays, such as chlorite, are modeled by randomly precipitating clays particles into the grain surface, already clay surface or quartz cement. The pore-filling clays tend to increase the tortuosity of the porous medium by clogging the pore throats.

Pore bridging clays, such as illite, are modeled as 3D lines that are near normal to the local direction of the skeleton [15]. These clays can affect the permeability of the medium building bridges throughout the pore throats and small pore bodies.

When feldspar dissolution occurs either partially or completely, this one is replaced by either void or clay.

The cementation process is the precipitation of minerals in the pore space. These minerals can be clays, carbonates or silica (quartz minerals). The role of the flow of pore fluid through the rocks is very important. If the flux is reasonably high, solids can be supplied in solution to the system and precipitation will cause a net porosity decrease without any loss in bulk volume.


Figure 2.2. 3D voxelized images, a) Pore space . b) Digital model of a sandstone rock, composed of quartz (yellow), clays (brown) and pore space (blue).

The end result of the the digital pore scale modeling can be observed in figure 2.2, a 3D voxelized images from a sandstone digital model.

2.1.2. Statistical Reconstruction

The microstructure characterizatrion is an important factor to obtain valuable predictive information from reservoir rocks [66]. Once the digital rock model is created, a set of statistical tools are available to describe the microstructure of the model in terms of spatial distribution of pore space and their connectivity. Two point Correlation function

The two point correlation function (equation 2.3) will find the possibility that in the same phase of material (pore or matrix), two points at different distance can be found [67].

The phase function Z(x) describes the medium indicating when a point x is inside of the pore space (Z(x)=1).

The medium can be described by the phase function as:

1 if x Ñ” pore space

Z(x) =

0 otherwise (matrix)

The intragranular porosity can be calculated as follows:


The point at which the correlation function is equal to zero is called the autocorrelation length or decay length. It does not mean that at this length the phase function is zero, it means that there is no spatial correlation anymore.


Figure 2.3. Two point correlation function of a digital model of a Fontainebleau sandstone rock. The decay length a is shown.

In the figure 2.3 the side length a or decay length is indicated, which represents the maximum grain size encountered in the model.

Lineal Path

The lineal path function PL(u) is defined in equation 2.4 and represents the likelihood of having a line segment of a given length and orientation in the phase of interest when randomly thrown into the sample. Once quantifying the l(u) segments that falls completely into the pore space, PL(u) is evaluated in every orthogonal direction for each pore voxel [14].

The parameter m denotes the quantity of placements of the segment of length u.

Figure 2.4. Lineal path function of a digital model of a Fontainebleau sandstone rock

Local Porosity Distribution

The local porosity distribution is defined in equation 2.6 [25].

A measurement cell φ(x, L) represents a cube centered at the vector x that has a side length L. A group measurement cells could obtain porosity as a local geometric properties in the following way:

Therefore, the local porosity distribution is obtained by

with the indicator

where µ(φ,L) is the empirical probability density function of local porosities.

Figure 2.5. Local porosity function of a digital model of a Bentheim sandstone rock

Local percolation function

The local percolation function characterizes or represents the connectivity of measurements cells of a given [25]. The connectivity of the pore space controls the transport in porous media.

The total fraction of cells, P(L) is given by

where Λ(x,L) is an indicator of percolation

Figure 2.6. Local percolation function a digital model of a Fontainebleau sandstone rock.

Statistical Procedure

The procedure implemented using the statistical tools in this work are as follow:

The two point correlation function is acquired and the decay length or autocorrelation length a is determined. Its value is an indication of the maximum grain size of the digital sandstone rock. In addition, if the two point correlation curve have more than one crossing point with the G(u) = 0, then it might be related with the presence of anisotropy within the microstructure of the digital sandstone rock.

The side length L of the measurement cell φ(x, L) or measuring box is set to be equal to the autocorrelation length a and m determine the smoothness of the local porosity distribution curve. This curve is computed and determined in a graph where the peak represents the most likely porosity value of the digital rock.

The same m value used in the local porosity distribution function is used in the lineal path function. At the point where the function PL(u)=0, the value of the x axis will determined the probability that a segment of line could fall entirely in a pore space of the size x.

Using the values of L and m previously defined, the local percolation function is obtained and it is represented in a graph. The local percolation curves for the x, y axis are found. The steeper the curves are ,the greater the connectivity of the pore microstructure is.

2.2. Pore Network Modeling

The following section reviews theories of author's such as Øren et al. [12, 13, 14, 15], Adler [26] , Hilditch [27], Blunt et al. [ 9, 31 ], Lake [30] , Vizika et al. [32].

Although is possible nowadays to compute single-phase and multiphase flow properties straight on the 3D voxelized image, the latter are computationally costly and they are limited to small 3D models. Therefore, network models are used as a simplification to simulate multiphase flow through idealized representations of the pore space. The models represent the microstructure of the pore space and flow processes. To describe the pore space in the network, a complete characterization of the 3-D structure is described by the pore structure, pore geometry and spatial continuity [12].

In the pore structure, the skeleton provides the necessary information regarding the pore space connectivity [12].

The thinning algorithm is the process used in order to extract the skeleton from binary 2D images and it is based on eliminating or preserving the interest pixels of the surrounding 8 neighbor pixels. In contrast, in 3-D, thinning algorithms used a collection of 0s and 1s combined and based on 26 neighbors per voxel. Another alternative is the inverse function of thinning called as dilation. Consequently, the ultimate dilation of the complementary grain network is precisely as the ultimate thinning of the pore network [12].

The 3-D skeleton process starts by first numbering each grain that formed the new numbered grainbed. Subsequently, this new grainbed will passed through an ultimate dilation of its grains.

The skeleton of a pore network is ultimately formed and fairly approximated by the ultimate dilation method [12] and its completion relay on the connectivity of its nodes.

The 3D voxelized image or geological model is represented by pore bodies and pore throats. Keeping the relevant features of the voxelized image structure for fluid flow in the pore network, the nodes and links are introduced. These represent a geometrical idealization of the pore bodies and pore throats respectively. The procedure is carried out using the following descriptors [12]:

(i) the maximum radius recorded inside a pore body defines a "node".

(ii) The smallest radius recorded in the pore space interconnecting two pore bodies inside a pore throat defines a "link".

(iii) The shape of a pore (defining multiphase flow in the same pore).

Once the descriptors are defined, the mapping procedure starts:

Each pore throat located within the voxelized image will be mapped using a plane normal to its wall. A vector rotating in increments of 10° 36 times is used to determined the radius or the detachment of skeleton voxel until the pore throat wall, see figure 2.7.


Figure 2.7. Principles of geometrical analysis of pore throats. a) skeleton line and the planes for the geometrical analysis. b) Cross sectional plane.

In the same way the geometry of the pore bodies are determined. In addition, the volume of the pore body is determined assigning an average radius from its centre point from the skeleton node [12].

An important parameter that defines the idealized cross-sectional shape of a pore body and throat is called shape factor, G. It is defined as G = A/p2 and its possible values ranges from 0.0481 (equilateral triangle) and zero (slit shaped). The shape factor is composed by the area A of the pore body (throat) cross section and the its perimeter lenght p2 [12].

The general idea is to find an idealized shape that substitutes the irregular shape of a pore body or pore throat of the digital rock using simple elements such as circular, square or triangular cross-section. The process starts determining first the shape factor G of a pore throat or pore body from the reconstructed pore space. Subsequently, the idealized irregular triangle cross-section is obtained using the same shape factor G, see figure 2.8.

Figure 2.8. Representation of a pore body (from voxelized image) with its respective shape factor G into an idealized figure. The idealized shape is an irregular triangular shape where β1, β2 and β3 are the half angle of each corner.

The equation used from elementary geometry that relates the half angle β of the triangle idealized shape with the shape fact G is the following:

In order to use equation 2.10, two conditions exist:

If G=0.0481 (equilateral triangle) then β1= β2= β3 and equation 2.10 can be verified.

If G < 0.0481 then from a single value of G different shapes can be obtained. Therefore, for a given G, β2 is randomly picked from a range set by β2,min ≤ β2 ≤ β2,max, where β2,min represents the case for an isosceles triangle (assuming β1= β2 ) and β2,max represents the case for a scalene triangle .


Using the random number β2 and substituting it in the equation 2.11, the values of β1 and β2 can be calculated. The other half angle can be calculated from equation 2.10 [14].

An important topological feature of a pore network is the so called coordination number, Z [28, 29]. The coordination number of a node is determined knowing the number of links connected to it and describes the connectivity of the pore network. According to P.E. Øren and S. Bakke, for a Berea sandstone digital rock, the coordination number ranges from 1 to 16 with an middling value of 4.45 [14]. Figure 2.9 shows a representation of the pore body and pore throat and the coordination number.

Figure 2.9. Pore body and pore throat using a "ball and stick" representation and their respective coordination number.

The importance of the coordination number relies on their influence on the flow behavior of the network and the importance of the residual non-wetting phase getting trapped [30, 31].

The nodes and links represent the pore bodies and the pore throats of a pore network. The effective node lengths and link lengths are defined by equation 2.12 and 2.13 respectively, where is the distance between the channel constriction and the node i (j), ri (j) is the radius of the pore body i (j) and rt is the radius of the constriction. The total length between two connecting nodes is defined by the length of the pore channel minus the length of the effective pore body and throat length.

Figure 2.10. Representation of two pore bodies i and j connecting by a pore link.

Thus the pore volume of the network is calculated by,

where ltot and ntot are the total number of pore throats l and pore bodies i respectively, Vp is the inter-granular porosity and Vt is the inter-granular porosity for pore and bodies respectively. Furthermore φµc represents the clay microporosity and Vc,i(j) are the clay counted in every pores and throats [13, 14].

Figure 2.11 Pore network model extracted from the geological pore scale model.

In order to model clay microporosity in the network, fractal cross sections with parallel pores are used. The model is described by Vizika and Lenormand [32] where an iterative process in which the perimeter of a circular pore ro is divided by 2 and subsequently in multiple parts and subsequently changing each part by half of the circle. The radius ro = ½.a.φµc where a is the resolution of the reconstructed model. The linear fractal dimension DL=1.6, used commonly for sandstone [14].

2.3. Network Fluid Modeling

The network fluid modeling aim is to simulate two-phase flow in the pore networks. It is assumed extremely low flow rate. Viscous forces are neglected across the network and therefore the capillary pressure is the dominant force within the fluid configuration.

The displacing fluid is injecting through every pore throat and every pore body on the inlet side, and the fluids are produced through the opposite side of the inlet face.

Periodical boundary conditions are imposed along the sides parallel to the main direction of the flow and it is assumed that the system is held at some constant reference pressure.

Primary Drainage

In the beginning, the network is 100% saturated with water (all the pores and throats), therefore it is strongly water-wet. Oil enters the network and reaches the available pore bodies and pore throats with the lowest threshold capillary pressure. This is forming the basis of the invasion percolation algorithm [29, 33]. Threshold or entry capillary pressures are calculated with the Mayer and Stowe [34] and Princen [35, 36] (MS-P) method. The threshold capillary pressure is governed by both pore shape and the receding contact angle and as well parameters as interfacial tension and radius using the following equation,

where r is the inscribed radius of the pore body or throat, θr is the receding contact angle, G is the pore shape factor, Fd is a function independent on the particular corner angles and therefore is not universal for a specific G. For strongly water-wet systems (i.e. ), then Fd = 1.

For threshold pressure, the presented method is generalized to the model introduced by Mason and Morrow [37, 38 ] and applied to any shape of triangular pores. The goal is to derive a general expression in terms of G (pore shape factor) for drainage. These computations are similar to the work done by Ma et al [39] for the displacement curvature of a meniscus in regular polygonal tubes.

At the threshold capillary pressure, oil enters to the throat with a fixed curvature and displaces water from the central part of the throat, leaving some of the water as AM's (Arc Meniscus) in the corners. In the absence of gravity the curvature of AM's is the same as the curvature of the invading interface.

Displacing a small distance dx of the AM's, then the work of the displacement must be balanced by the surface free energy change, see equation 2.17.

where the subscript s represents the solid, Aeff the effective area occupied by oil, LOS the length of the solid wall in contact with the oil and LOW the perimeter length of the arc meniscus.

Using Young's equation , where , and are the interfacial tension of the oil-solid , water-solid and oil-water interface respectively with equation 2.17, equation 2.18 is obtained [13].

Therefore, in the eq. 2.21 the curvature of the invading interface at the condition for invasion is given as the ratio Leff/Aeff,, Parameters Aeff, LOS and LOW are determined from elementary geometry and are obtained by equations 2.19, 2.20 and 2.21 respectively. The following equations for primary drainage are taken from the work developed by Øren et al. [13].

where rd is the radius of curvature, is the corner half angle of the triangular cross-section and S is the perimeter.

The following equation is obtained after some algebraic operation in order to obtain a quadratic expression for rd that has a solution:

where . The equation 2.22, could expressed as a function of the ratio r/rd

where the function is given by,

If AM's is present in all the corners, the expression for D becomes,

In this way D is dependant only on θr and .

The purpose of the series of equations shown from 2.19 until 2.23 is to derive a function (eq. 2.24) that relates the receding contact angle with the curve shape factor G. This function will play an important role in the drainage capillary pressure threshold.

Primary Imbibition

The oil pressure is kept it constant at the outlet whereas the water is slowly being injected into the system. The capillary pressure drops during water injection.

Three types of displacement mechanisms were described by Lenormand et al.[6], namely piston-like, pore body filling and snap-off. Below, main characteristics are retracted from the theory of Øren.

During the piston type displacement, the oil is displaced from a pore throat by an invading interface initially located in the nearest water filled pore body.

The oil and water interfaces of the arc meniscus in the throat remain pinned at the position Lb that was established at the end of primary drainage.

where Lbi is the distance between the pinned contact line and the corner i and . The hinging angle and the relationship with and rpt can be determined by the following equation:

The effective perimeter Leff is given by

where the angle is given by

The equations 2.28, 2.29 and 2.30 constitute a non linear system that can be solved numerically for the radius of curvature rpt and the threshold capillary pressure .

Solution for the normalized threshold capillary pressure shows that spontaneous imbibition by piston type displacement always occur for , less than 90° and may occur for greater than 90° [13].

Using the equation 2.29, the maximum advancing angle for this process is acquired [13]

If the advancing angle is greater than , the water invasion is forced. Forced water invasion is analogous to the oil drainage process. If , the entry capillary pressure may be determined from Eq. 2.29. If , the threshold capillary pressure is given by

The pore body filling is a displacement mechanism that depends on the pore body dimension and the pore throat collection that are filled with oil [13].

In the case that a pore body has a coordination number z, there are z -1 such mechanisms, referred to as I1 to Iz-1. If on throat that is connected carry oil (i.e., I1), then the displacement method is equivalent as the piston type of displacement and so the threshold capillary pressure. If there is more than one connecting throat, then more complex mechanism are formed.

If t he mean radius of curvature for filling by an In mechanism is computed by

where rp is the pore body radius , bi is the input parameters, ri is the radii of the oil filled throats, and xi is the random numbers between zero and one. The water invasion is forced when the I1 mechanism or piston type displacement occurs and the capillary pressure threshold is determined by . If [13].

The last displacement method used is known as Snap-off. [13].

For a strongly water-wet system, the threshold capillary pressure is given by [13].

Once the capillary pressure decreases to certain levels and there is the presence of the contact angle hysteresis, the AM's advance to the center of the pore space.

If , snap-off displacements occurs at a positive capillary pressure[13]. The capillary pressure at which this occurs can be calculated using the following equation:

When spontaneous imbibition occurs by snap-off. If, the invasion is forced. The capillary pressure threshold for this process it is given by

In case of larger values of , the AM's in the corner remain pinned whereas the two others advance at . The snap-off occurs when either when the two advancing AM's meet or when the AM's in the sharpest cornet meet the pinned AM. The threshold capillary pressure for the latter case is given by,

Piston like advance in a mixed-wet pore body or throat leaves the center of the pore space filled with water. If, a film of oil may be left sandwiched between the water in the center and the water in the corner(s) , refer to figure 2.12. Therefore the formation of oil films falls in the increase of the connectivity of the oil phase.

Figure 2.12. Water occupying the center and corner of a pore body or throat whereas the oil is sandwiched between them. Triangular pore section.

Once the oil/water interface meet on the film, the stability will be lost and the equation for the capillary pressure sould be the following:


In this work, wettability effects at pore-scale level are shown in two different processes: primary drainage and first imbibition.

Primary Drainage

Pores and throats are considered with different cross-sectional shapes. The filling procedure is done in order of capillary pressure. In figure 2.13 can be seen the effect when the water remains at the corner and the oil in the center of the triangular pore.

Figure 2.13. Water and oil in a pore body or throat with triangular cross-section. After primary drainage oil enters in the pore space while the water remains in the corner.

Primary Imbibition

At this stage, the oil has already invaded the pore body or throat that was completely filled by water. The water forms a stable film that covers the pore structure until a critical capillary pressure is reached. Once it has occurred, the protective film collapse and the wettability of the pore structure changes [40]. According to Kovscek et al. [40] the critical capillary pressure for a pore body or throat i, is giving by

Where П represents the film curvature, Г is the disjoining pressure. Both are input parameters and x varies uniformly between 0 and 1. Primary drainage terminate until some specific , or minimum water saturation is reached.

The water saturation can never be zero. The reason for this is that clays are modeled having micro-porosity which remain water saturated during the displacement process [14].

During this section, the analytical model to account for mixed wettability will be explained. It was proposed first by Dixit et al. [41] deriving an analytical relationship between IAH and α1 (fraction of oil-invaded pores) for a group of capillary tubes. In a later stage Øren and Bakke[14] extended their analysis. They included effects such as: accessibility effects, contact angle and phase trapping.

The very important parameters in this model [14] are α1 that oil invades pores that become oil wet (oil-wet pore fraction), contact angles and the spatial distribution of oil wet pores and.

The Amott water index, can be expressed as:

where ΔS1 is the water saturation change due to spontaneous imbibition. ΔS2 is known as the extra water saturation change accounting for forced displacement. In the same way, the Amott oil index can be expressed as,

For the first spontaneous water imbibition, the saturation change can be expressed as :

where Sow1 is the oil wet pores volume fraction, occurred when (θa>90°) and a1 is the water wet region trapping function defined as follows,

where is the residual oil saturation in water pores, is the volume fraction of oil wet pores and defined as , depends on oil wet pores distribution. According to Øren and Bakke [14] three pore size distributions are considered:

Randomly distributed oil wet pore bodies (MWR)

Large pore bodies preferentially becomes oil wet (MWL)

Small pore bodies preferentially becomes oil wet (MWS)

In case the oil-wet pore are randomly distributed, then . In contrast, the saturation function,, varies almost linearly with such that if the oil-wet pores are not randomly distributed. This behavior was obtained by Øren work [14].

Conversely, for the force displacement, the change in water saturation is found with the following equation:

where a2 is the same type of function as a1 but for the oil region.

where is the residual oil saturation in the oil wet pores. Using equations 2.40, 2.42 and substituting in

equation 2.38, the fraction of oil wet pore can be obtained in terms of Iw, as

secondary drainage is out of the scope of this work, therefore for more details reader can refer to[14].

It was demonstrated by Øren that when α1 > α1c (Iw<1), the value of the trapping function a1 is fairly close to 1.

An approximate expression for a2 was found when plotted (1-a2) versus Sow1c/Sow1 where a power law behavior was observed and a following empirical expression was established,

The contact angle hysteresis arises from the difference between the effective advancing contact angle during waterflooding and the effective receding contact angle during drainage. This results from the roughness of the pore surface, pore geometry effects and absorption [42]. To account for the hysteresis, and are characterized in terms of the equilibrium angle where and are defined as equilibrium angles at and respectively. In this model it is assumed that , for the water-wet pores, to be randomly distributed . Correspondingly, for oil-wet pores, the contact angle distribution range is defined from and that depends on the fraction of oil wet pores during secondary drainage α2 and given by

Transport properties

Absolute Permeability

Predicting the permeability of real porous media has obvious practical significance. However, so far it has been difficult to determine this parameter.

Nevertheless, if the medium is properly characterized then flow assumptions can be made. Fluids are considered to be Newtonian, incompressible, immiscible, presenting laminar flow and thereafter for a low Reynolds number the flow of the fluids involved are governed by Stokes equation:

where v, p and µ are the velocity, pressure, and fluid viscosity , respectively. In general v satisifies the conditions. The equation 2.47 arises from mass conservation and the equation 2.48 derives from the momentum.


is partially periodic

and Sp represents the surface of the wetted solid inside the unit cell .

In order to solve previous Stokes equations, a D3Q19 Lattice Bolztmann algorithm [43] is used.

For this study, the method used is known as relaxation Boltzmann scheme or 'BGK' (Bhatnagar-Gross-Krook) method [43]. In this method, at each time step the particle distributions propagate from one lattice node to the neighboring ones redistributing locally in subject to the conservation of mass and momentum by a collision operator. The mathematical expression for collision and propagation is shows as follows:

where Fi(x, t) considered as the probability of finding a particle with the lattice velocity ci, in the location x, and at time t. The relaxation parameter ω determines transport properties of the fluid.

To simulate fluid flow, the lattice is identical as the discretized rock image or the voxel image. The example is shown in figure 2.14. Each lattice node is put in the center of a voxel in the pore space image.

The hydrodynamic fields, such as the local macroscopic mass density ρ and the fluid velocity v, are defined for a given node in terms of moments of the particle distributions Fi(x, t) by

where the number of admissible lattice velocities per node is denoted by n. It is referred to the particles at each node that can move along n different directions or stay still. The number of admissible velocities at each node is related to the connectivity assumption adopted in the simulations. The standard connectivity numbers are 6, 18 (where two voxels are linked if they have a common face or a common edge) or 26.

Figure 2.14. Lattice structure of the lattice Boltzmann model with 19 velocity directions per lattice nodes

As described in fig. 2.14, each node is coupled to its six nearest and twelve diagonal neighbors at a distance of , where a is the voxel resolution. Particles can only stay at the node or move along the fixed links connecting the lattice nodes, and interact at the nodes. All particles associated with a given direction have the same lattice velocity.

In order to calculate the absolute permeability of the porous medium in steady state flow, the initial conditions in the simulation are set v = 0 for the flow velocity and ρ = 1 for the fluid density in the whole domain. The LBM use a no-slip boundary condition at the solid-fluid interfaces, and a prescribed flow at the inlet and outlet faces. To simulate the no-slip boundary condition at the solid-fluid interfaces, a complete bounce-back scheme is used [43, 44, 45].

Here, a constant and consistent body force in order to drive the flow is used. The uniform body force is intended to create a momentum into the flow as the true pressure gradient ΔP/L1. To implement this scheme, a constant value ΔFi(x, t) was added to the particle distributions moving along the pressure gradient, and place a corresponding penalty on the counter-particle distributions

Assuming that the magnitude of ΔFi(x, t) be equal for all these directions, we have

In addition, a periodic boundary condition is applied at the inlet and outlet to allow the departing particle distributions to reenter the flow domain through the inlet.

The single-phase, low Reynolds number flow through a porous rock sample is described by Darcy's law, through the following equation:

where k_ is the permeability of the porous medium, is the fluid density, is the fluid kinematic viscosity, ΔP/L1 is the total pressure drop along the sample length L1, and q is the volumetric fluid flux through the porous medium,. Here φ is the sample porosity, and is the velocity averaged over all lattice nodes in the pore space. Using the scaling forms mentioned above, and Eqs. (2.53), (2.52) and (2.51), the permeability of the sample can be calculated directly from the following equations,

In a LB simulation, only in the connected (percolating) part of the pore space, fluid can flow through the sample. In other words, the disconnected pore space does not contribute to the sample permeability. Identification of the percolating pore space is achieved by clustering the pore space using the Hoshen-Kopelman algorithm.

Once the permeability k is scaled back an average of the x,y,z directions are calculated .

The procedure to obtain absolute permeability was obtained from different works [12, 43]

Formation Factor

The formation factor F is defined as

where is the conductivity of the porous medium and the bulk electrical conductivity of the fluid that fills the pore space.

The formation factor is computed independently in all three axes. As can be seen from equation 2.58, there is a linear relation between the total potential drop ΔФ over the distance Δx in the direction of the flow and the total electrical flux E.

where A is there total cross-sectional area normal to the flow direction. Potential gradient is applied across the sample and no flow boundary conditions are assumed on the others. The electrical potential is determined by solving the Laplace equation:

where n is the unit vector normal to the pore wall Sp. For sandstones that precisely do not contain clay, the assumption of an insulating phase is verified. The microscopic fluid conductivity is given by

The flux E at any cross-section

for any x in the domain. From equation 2.58 and 2.64 can be obtained the following expression for the electrical conductivity.

and then formation factor .

Equations 2.57 and 2.60 are solved by a finite volume method on the regular grid that follows from the binary representation of the sample. Equation 3.76 then takes the form

If a neighbouring control volume is closed (i.e. matrix) then the boundary condition in equation 2.60 determined that there is no contribution from the face in equation 2.63. Due to the implemented Dirichlet boundary condition by

The resulting Ax = b system can be resolved by the Conjugate Gradient method equipped with and SSOR (Symmetric Successive Over Relaxation) preconditioner [46].

The method works as follows:

The symmetric matrix A is decomposed as follows:

where L is the strictly lower part of A and D is the diagonal of A.

SSOR Preconditioner [46] consists in taking the equation 2.68 and used the relaxation parameter fixed at the interval [0, 2] as a sufficient condition of this SSOR method algorithm.

In such a way, the unknowns are calculated in order to obtain and therefore compute Formation factor in each direction.

Calculating the harmonic mean of the formation factor directional dependant will give the average formation factor and is determined as follows:

If the medium is isotropic Fx = Fy = Fz = F. Otherwise, the degree of anisotropy, is expressed in terms of heterogeneity factor Af which is defined as

The procedure implemented to obtain the formation factor was obtained from different works [12, 13, 47].

2.5. Constitutive relations

First, fluid areas and conductance has to be calculated. A triangular cross-section area is assumed in pore bodies and pore throats and varies depending if the arc meniscus in the corners of pore bodies and pore throats come from water or oil. In case the area is occupied by water then,

If is the oil is present in the corner of pore bodies and pore throats as a film then,,

Assuming Poiseuille law, the hydraulic conductance of completely filled pore body and throats is computed in the following way;

where q is the flow rate, g is the hydraulic conductance and is the pressure gradient. Assuming incompressible flow and boundary conditions along the walls of a pore body or a pore throat, hydraulic condutances are determined by integrating the velocity field over the pore body or pore throat cross section [13]. Therefore, the following equation is determined,

Calculating between two connecting pore bodies (nodes) I and J the flow rate. (Laminar flow is assumed).

where LIJ is the spacing between the center of the nodes, Δpi,IJ is the pressure drop over the element, qi,IJ is the flow rate and gi,IJ is the hydraulic conductance of fluid i between the centers of the nodes I and J. The hydraulic conductance gi,IJ is considered to be calculated as the harmonic mean between the link and the two-half nodes and is calculated in eq. 2.76.

where t is the subscript for pore throats.

Applying mass conservation for each node i:

where J pass through all the pore thoats connected to the Pore Body I. Using equation 2.72 and 2.74 will give a set of linear algebraic equations of node center pressures. Applying a constant pressure is used in the inlet and outlet. The expression can be written as a matrix,

where G is a sparse matrix whose element's contains the conductances, P is the nodal pressure(unknown) and b is a vector containing zeros except for the inlet and outlet and interface nodes.

In order to calculate the permeability of the network, a constant pressure gradient along the system is applied. The eq. 2.78 term P, give the values of the nodal pressures that can be solved using the conjugate gradient method [46]. Finally, from P, the flow rate is calculates and thus the absolute permeability using Darcy's law.

Relative permeabilities are calculated from Darcy's law as well. The pressure drop is calculated in each phase during the different stages of the displacement in which saturation and total flow rate are determined [13].

Due to the analogy between Poseuille's law and Ohm's law, the flow of electrical current in the network is also described using eq. 2.73 and 2.75 where the pressure is replaced by voltage and hydraulic conductance by electrical conductance.

The electrical conductance ge assumptions:

1) It is geometry dependable, i.e. pore body or throat.

2) The walls of the pore space are insulated.

Putting the assumptions in place the following equation is shown,

where is the electrical conductivity of the fluid.

The electrical conductivity, σ, of the pore network is computed by applying a constant voltage drop across the network and then applying eq. 2.74 to find the flow of current between each pore body. Finally, the total current is calculate in order to determine the , see eq. 2.80.