**Percolation Approach in Underground Reservoir Modeling**

Mohsen Masihi1 and Peter R. King2 *1Department of Chemical and Petroleum Engineering, Sharif University of Technology, 2Earth Science and Engineering Department, Imperial College London, 1Iran 2UK* 

**1. Introduction** 

262 Water Resources Management and Modeling

Choubey, V.K. (1996) Assessment of waterlogged area in IGNP Stage-I by remotely sensed

Dafny, E.; Burg, A. and Gvirtzman, H. (2010) Effects of Karst and geological structure on

Ehlers, M.; Edwards, G. and Bedard, Y. (1989) Integration of Remote sensing with

Indo-Dutch Network Project-IDNP (2002) A Methodology for Identification of Waterlogging

Lillesand, T.M. and Kiefer, R.W. (2004) Remote Sensing and Image Interpretation, 5th

Middlemis (2000) Groundwater Flow Modeling Guidelines, Murray-Darling basin commission. Aquaterra Consulting Pvt, Ltd, S. Perth, Western Australia Peleg, N. and Gvirtzman, H. (2010) Groundwater flow modeling of two-levels perched

PPSGDP Report (2000) Groundwater Flow and Solute Transport Numerical Models For

Qureshi, M.K.A.; Sheikh, M.I. and Khan, A. (2003) A Report on the Geology and Environmental

Sarwar, A. (1999) Development of a conjunctive use model, an integrated approach of

Shelton, R.L. and Estes, J.E. (1980) Remote Sensing and Geographic Information Systems: An

Siddiqi, J. (1999) Effect of El-Nino on Agriculture of Pakistan, Pakistan Agricultural

Soil Survey Report (1967) Reconnaissance soil survey of Gujrat district, Directorate of Soil

Thangarajan, M. (2004) Regional Groundwater Modeling, Capital Publishing Company,

UNESCO (1999) Water Resources of Hard rock Aquifers in Arid and Semi-arid Zone. Ed. J

WAPDA (1993) Pakistan Drainage Sector Environmental Assessment – National Drainage

Yang, J.W. and Radulescu, M. (2006) Paleo-fluid flow and heat transport at 1575 Ma over an

Zaisheng, H.; Wang Hao and Chai Rui (2006) Transboundary Aquifers in Asia with Special

Programme Vol. 4: (Data (water, Soil and Agriculture), NESPAK and Mot

E–W section in the Northern Lawn Hill Platform, Australia: Theoretical results from finite element modeling. Journal of Geochemical Exploration, 89(1-3): 445–449.

groundwater flow: The case of Yarqon-Taninim Aquifer, Israel. Journal of

Geographic Information Systems: A Necessary Evolution. Photogrammetric

and Soil Salinity Conditions Using Remote sesning. CSSRI, Kamal and Alterra-

karstic leaking aquifers as a tool for estimating recharge and hydraulic parameters.

SCARP II and Shahpur Area of Chaj Doab, Tech. Rep. No. 39, Punjab Private Sector Groundwater Development Project Consultants, Project Management Unit,

concerns of the Rasul-Mandi Bahaudin Quadrangle, Punjab, Pakistan, Ministry of Petroleum and Natural Resources, Geological Survey of Pakistan, Quetta, Pakistan. Russo, S.L. and Civita, M.V. (2009) Open-loop groundwater heat pumps development for

surface and Groundwater modeling using a Geographic Information System (GIS).

and field data. Hydrology Journal, Vol. XIX (2), pp. 81-93.

Doherty, J. (1995) PEST ver. 2.04. Water mark Computing, Corinda, Australia

Engineering and Remote sensing, 55, no. 11: 1619-1627.

Irrigation and Power Department, Govt. of Punjab, Lahore

large buildings: A case study. Geothermics, 38(3): 335–345.

Unrealized Potential. Geo-processing, 1, no. 4: 395-420.

MacDonald Int. Ltd., WAPDA, Lahore, Pakistan.

Emphasis to China, UNESCO Report

Hydrology, 389(3–4): 260-275.

ILRI, Wageningen: pp. 78

Edition, John Wiley & Sons.

Journal of Hydrology, 388(1-2): 13–27.

PhD thesis, University of Bonn, Germany.

Research Council, Islamabad.

Survey of Pakistan, Lahore.

New Delhi, India.

W Lloyd.

The spatial distribution of the underground heterogeneities which may be appeared on various scales can affect the flow and transport of fluids (e.g. underground spread of pollution or recovery from hydrocarbon reservoirs). Reservoir modeling is conventionally used to incorporate the effect of uncertainty in the spatial distribution of geological heterogeneities (such as fractures as appeared on various scales as emphasized by Odling et al, 1999) and to assess the reservoir performance parameters such as the reservoir connectivity and conductivity. It needs detail data to make geological model, then upscale it to coarser grid and run flow simulation. Since the natural variability of geological formations extends over many length scales, this approach requires very fine modeling schemes. In addition, the predictions may not be achievable in the early life of reservoirs because of the sparse data. Moreover, in stochastic modeling in order to make a reliable underground prediction or for sensitivity purposes, it is necessary to construct a number of possible reservoir models. The problem with this approach is that it is computationally very expensive and time consuming. Hence, there is a great incentive to produce much simpler physically-based model to predict the underground reservoir performance very quickly.

In this chapter, we used the idea that the flow characteristics in porous media are mainly controlled by the continuity of conductivity contrasts (e.g. high hydraulic conductive streaks) to model stochastically the underground reservoir flow. There are many cases for both hydrologist and petroleum engineers where the geological formations consists of a mixture of good sandstones with high hydraulic conductivity (i.e. flow units) and poorer rock types (e.g. siltstones, mudstones and shales) with negligible permeability. Good sandstones with a higher permeability and porosity are the main bodies containing fluid within their pores. Obviously the volume of fluid in place and recoverable fluid between a pair of wells depends on the connected fraction of these good fluid bearing zones. Examples include fluvial sediments containing paleo-channels (high conductivity zones embedded in a low conductive background), shale/sandstone sequences (non-conductive inclusions embedded in a conductive matrix), fractured formations (with fractures as the connected high conductive zones), and coastal deposits (deltaic systems representing the conductive media).

As the mathematical model of connectivity and conductivity in complex geometries is given by percolation theory (Stauffer and Aharony 1992) we use "percolation approach in reservoir modeling" for the title of the presented method. Imagine a typical reservoir model to be constructed with an object based technique. Assume a simple permeability model of the reservoir (black and white model for sand/shale or fracture/matrix). Then consider sandbodies/fractures as simple geometrical objects located in space with simple statistics and use percolation theory to estimate the uncertainty in reservoir performance. In particular, we can estimate the probability that wells are in connection via the sands, the fraction of the total sands which is in contact with these wells, the reservoir conductivity, potential oil recovery, sweep efficiency, breakthrough time and post breakthrough behavior within a particular well configuration. Throughout the field life many decisions have to be made based on the technical factors (see Table 1) which can be addressed by using percolation approach.


Table 1. A summary of the key decisions and technical factors during underground hydrocarbon reservoir life

Moreover, the percolation approach is able to estimate the uncertainty in reservoir parameters which is not possible with a single realization reservoir model. The advantage is that the effects of the complex geometry which influence the flow can be easily estimated in a fraction of a second on a spreadsheet by simple algebraic transformations (Masihi et al. 2007). Clearly the disadvantage is that much of the flow physics and subtleties of the heterogeneity distribution are missed.

### **2. Percolation theory**

The mathematical analysis of percolation theory was first developed by Broadbent & Hammersley in 1957 and since then the topic has been intensively studied [Stauffer & Aharony, 1992]. It has been applied in many fields from the spread of forest fires and the spread of diseases to the flow in porous media and fractured rocks [Sahimi 1994]. This theory links the global physical properties such as connectivity and conductivity to the number density of objects placed randomly in space and reduces many results to simple power laws from which all outcomes are predicted by simple transformations. A simple percolation model is an infinite size lattice of squared sites (so called *site* percolation) that are occupied randomly and independently to an occupancy fraction *"p".* Neighboring

As the mathematical model of connectivity and conductivity in complex geometries is given by percolation theory (Stauffer and Aharony 1992) we use "percolation approach in reservoir modeling" for the title of the presented method. Imagine a typical reservoir model to be constructed with an object based technique. Assume a simple permeability model of the reservoir (black and white model for sand/shale or fracture/matrix). Then consider sandbodies/fractures as simple geometrical objects located in space with simple statistics and use percolation theory to estimate the uncertainty in reservoir performance. In particular, we can estimate the probability that wells are in connection via the sands, the fraction of the total sands which is in contact with these wells, the reservoir conductivity, potential oil recovery, sweep efficiency, breakthrough time and post breakthrough behavior within a particular well configuration. Throughout the field life many decisions have to be made based on the technical factors (see Table 1) which can be addressed by using

Phase Business Decision Technical Factors

Plateau End of plateau Breakthrough time and post

Table 1. A summary of the key decisions and technical factors during underground

Moreover, the percolation approach is able to estimate the uncertainty in reservoir parameters which is not possible with a single realization reservoir model. The advantage is that the effects of the complex geometry which influence the flow can be easily estimated in a fraction of a second on a spreadsheet by simple algebraic transformations (Masihi et al. 2007). Clearly the disadvantage is that much of the flow physics and subtleties of the

The mathematical analysis of percolation theory was first developed by Broadbent & Hammersley in 1957 and since then the topic has been intensively studied [Stauffer & Aharony, 1992]. It has been applied in many fields from the spread of forest fires and the spread of diseases to the flow in porous media and fractured rocks [Sahimi 1994]. This theory links the global physical properties such as connectivity and conductivity to the number density of objects placed randomly in space and reduces many results to simple power laws from which all outcomes are predicted by simple transformations. A simple percolation model is an infinite size lattice of squared sites (so called *site* percolation) that are occupied randomly and independently to an occupancy fraction *"p".* Neighboring

Connected oil volume fractions

Flow path tortuosity and connected

Distribution of remaining oil pools

connected by wells Swept volumes Effective permeability

breakthrough behavior

volume of oil

Number and location of wells

Reserves calculation Production rates

Decline Decline rate or increase of water or gas Infill locations

percolation approach.

hydrocarbon reservoir life

**2. Percolation theory** 

heterogeneity distribution are missed.

Production

Exploration/ appraisal

sites are in the same cluster if they are both occupied (Fig. 1). This can be checked by an algorithm proposed by Hoshen and Kopelman 1976. Consider that sites are filled with a fluid then the cluster is a region accessible to the same well. Percolation theory describes mathematically how the number and the size of these clusters vary as a function of the "*p".* As this occupancy probability is increased, the finite size clusters bridged and grow in size. At one particular value (called percolation threshold,��) one large cluster (so called the spanning, percolating or infinite cluster) spans the whole region (the blue clusters in figures below). However, there exist also other small clusters which get absorbed as *p* further increases leaving a few isolated sites. From flow studies point of view across the region (as between two wells) we are interested in the spanning cluster. However, to analyse connectivity to a point (or a well) the contribution of the finite clusters then must be considered.

Fig. 1. Illustration of (left) various clusters in different colors before the threshold, and (right) spanning cluster in blue at the threshold

Other extensions of this simple site percolation model include using other lattice shapes (e.g. triangular, hexagonal or cubic lattices), other connection criteria (e.g. next nearest neighbors), directed percolation (e.g. contact is only in left-right direction), bond percolation, continuum percolation (i.e. geometrical objects are distributed randomly in space and correlated percolation (i.e. objects may not be distributed randomly). These do not alter the essential ideas of percolation theory; however, may alter some numerical values of the parameters used to describe the percolation. One of these adjusting parameters is percolation threshold which is dependent on the type/detail of the system. The percolation threshold �� *<sup>∞</sup>* value depends on the details of the system. The threshold value for some geometries in two and three dimensions are given in Table 1 [Stauffer and Aharony, 1992; Berkowitz and Balberg, 1993, Baker et al 2002]. As can be seen in Table 2 the threshold values in three dimensions are much lower than the values in two dimensional lattices. One fundamental percolation property is the connectivity �(�) which is defined as a probability that any point in space belongs to the spanning cluster which shows the strength of the percolating cluster and scaled as [Stauffer and Aharony, 1992]:

$$P(p) \ll (p - p\_c^{\infty})^{\beta} \tag{1}$$

where p� <sup>∞</sup> is the percolation threshold of an infinite system and the β is called the connectivity exponent. Above the threshold, the connectivity, �(�), increases rapidly. The other property is the correlation length �(�) defined as �� = ∑ ���(�)� ∑ �(�) (with �(�)as the two point correlation function or the probability that two sites separated by a distance � are in the same cluster) which represent a typical size of finite clusters and used in analysis of connectivity to a point (or well) [Stauffer&Aharony, 1992],

$$
\xi(p) \propto |p - p\_c^{\circ}|^{-v} \tag{2}
$$

With � as the correlation length exponent. Moreover, the mean cluster size��(�) follows a simple power law as,

$$S(p) \propto |p - p\_c^{\alpha}|^{-\nu} \tag{3}$$

Where the � is called the mean cluster size exponent.


Table 2. Threshold values of different shapes for site, bond or continuum percolation

Next property is the effective permeability. From percolation theory the effective permeability of an infinite system �(�)has a power law behavior [Stauffer and Aharony, 1992]:

$$K(p) \ll (p - p\_c^{\alpha})^{\mu} \tag{4}$$

where the universal exponent ��is called the conductivity exponent. Above the threshold �� *∞*, the effective permeability increases slowly as compared with the connectivity emphasizing that a part of the spanning cluster (called the backbone) contribute to the flow (Fig. 2). The backbone consists of blobs (i.e. the multiple connected part of the backbone) and the red links i.e. the link between two blobs (Fig. 2).

Fig. 2. Illustration of (left) the tortuous conducting path with many dead ends and (right) red links and blobs between two points on the backbone

the two point correlation function or the probability that two sites separated by a distance � are in the same cluster) which represent a typical size of finite clusters and used in analysis

�(�) ∝ |����

�(�) ∝ |����

Square 0.593 0.500 Triangular 0.500 0.347 Hexagonal 0.696 0.653 Simple Cubic 0.312 0.249

Table 2. Threshold values of different shapes for site, bond or continuum percolation

of an infinite system �(�)has a power law behavior [Stauffer and Aharony, 1992]:

With � as the correlation length exponent. Moreover, the mean cluster size��(�) follows a

�| ��

�| ��

Shape Site Bond Continuum

Disks 0.676 Aligned Squares 0.666 Spheres 0.289 Aligned Cubes 0.277

Next property is the effective permeability. From percolation theory the effective permeability

�(�) ∝ (����

 Fig. 2. Illustration of (left) the tortuous conducting path with many dead ends and (right)

where the universal exponent ��is called the conductivity exponent. Above the threshold ��

the effective permeability increases slowly as compared with the connectivity emphasizing that a part of the spanning cluster (called the backbone) contribute to the flow (Fig. 2). The backbone consists of blobs (i.e. the multiple connected part of the backbone) and the red

(2)

(3)

�)� (4)

*∞*,

of connectivity to a point (or well) [Stauffer&Aharony, 1992],

Where the � is called the mean cluster size exponent.

links i.e. the link between two blobs (Fig. 2).

red links and blobs between two points on the backbone

simple power law as,

Backbone between points *i* and *j* are the bonds that belong to at least one self-avoiding walk between *i* and *j*. The burning algorithm described by Stanley et al 1984 is used to identify the backbone. The backbone fraction of the percolating cluster is the amount of oil that can be swept during water/gas injection. The critical exponents �� �� �� � in these power laws are universal. By universality, we mean that their values are independent of the small scale structure (i.e. the kind of lattice, site or bond percolation, lattice or continuum systems). They only depend on the dimensionality of space (i.e. two or three dimensional systems). The generally accepted literature values for these exponents for two and three dimensional systems are given in Table 3 [Stauffer and Aharony, 1992; Berkowitz and Balberg, 1993]. The principle of universality is a very powerful concept as it enables us to study and understand the behaviour of a very wide range of systems without needing to worry so much about the fine scale details.


Table 3. Selected percolation exponents in two and three dimensions.

The basics percolation has been developed on lattices; while, in many natural systems the objects are not restricted to fixed points or the objects can be of any shape with variable length, direction and number of interconnected bonds. In such cases we use continuum percolation, where the geometrical objects are distributed independently and randomly in a region. Examples include fracture systems where there is theoretically no end to the degree of fracturing [Sahimi, 1995] or overlapping sandbodies with various shapes and sizes [Masihi and King 2008]. In continuum percolation models, the principal of universality allows us to use the same scaling laws with the same numerical values of the critical exponents as in the case of lattice percolation. However, we need to the percolation threshold as it depends on the details of the system. Examples of applying percolation theory to uncorrelated (or even correlated) continuum systems that check the universality and determination of the percolation threshold of different models can be found elsewhere [Gawlinski and Stanley, 1981; Lee and Torquato., 1990; King, 1990; Berkowitz, 1995; Lorenz and Ziff, 2001; Baker *et al*., 2002].

It should be emphasized that the percolation cluster and its other sub-networks have fractal properties. Mass of the percolating cluster (as defined as total number of the sites �) scale as, ����� where the fractal dimension �� of the percolating cluster is given by �� � � � �� � with *d* as the space dimension. A sub-network of percolating cluster is the Backbone which is a fractal object with its mass �� scaled with total number of the sites � as ������ (with �� ������ � where �� is the fractal dimension of the Backbone cluster and �� is the connectivity exponent of the backbone cluster). Another sub-network is called the red bonds which is also a fractal object with its mass ���� scaled with total number of the sites � as ���������� (with ���� � ��� where ���� is the fractal dimension of the red bonds). The numerical values of the fractal dimensions of various percolation networks are given in Table 4 [Stauffer and Aharony, 1992; Sahimi, 1995].


Table 4. The numerical values of the fractal dimensions of various percolation networks

The fractal behavior of percolation-based reservoirs can be observed during well testing results analysis. Conventional build up test for radial flow in a single well homogenous reservoir gives the pressure change �� over the time period��� as �����(��). However, for some field data the conventional model matches only a part of the data not being improved by considering other effects including formation damage along the fracture face or nondarcy flow near the wellbore. In such cases, the alternative model which gives a reasonable match is to use solution of pressure transient equations with rock properties have fractal behavior (Chang & Yortsos 1990).

#### **3. Two basic models**

Two main applications of percolation approach include low-to-intermediate net to gross reservoirs and fractured reservoirs. Geologists have used the idea that a reservoir consists of the geometrically complex connected and disconnected sandbodies to be able to qualitatively describe it by means of outcrop studies and subsurface and for the reliable estimation of accessible hydrocarbons in place and the expected recoverable oil [Haldorsen et al., 1988]. Hence, sandbodies within a reservoir are assumed to be flow units connected in a complicated way through the sedimentary processes along with and poorer siltstones, mudstones and shales (of lower permeability) [Bridge and Leeder, 1979]. An example is a meandering river which deposits sand layers over the time in its flowing bed results in formation of a system of embedded sandbodies in an impermeable background (Fig. 3). Moreover, outcrop studies often show a network of connected fractures (Fig. 3). In this study, we use overlapping sandbodies and fracture network models within the framework of continuum percolation to determine the reservoir connectivity and conductivity.

Fig. 3. Illustration of (left) Overlapping sandbodies deposited by a meandering river over the ages (Nurafza et al 2006). (right) 2-D fracture model (Odling et al 1999)

#### **3.1 Overlapping sandbody model**

268 Water Resources Management and Modeling

numerical values of the fractal dimensions of various percolation networks are given in

91/48 1.64 3/4

2.52 1.8 1.14

Fractal dimension 2D 3D

Table 4. The numerical values of the fractal dimensions of various percolation networks

The fractal behavior of percolation-based reservoirs can be observed during well testing results analysis. Conventional build up test for radial flow in a single well homogenous reservoir gives the pressure change �� over the time period��� as �����(��). However, for some field data the conventional model matches only a part of the data not being improved by considering other effects including formation damage along the fracture face or nondarcy flow near the wellbore. In such cases, the alternative model which gives a reasonable match is to use solution of pressure transient equations with rock properties have fractal behavior

Two main applications of percolation approach include low-to-intermediate net to gross reservoirs and fractured reservoirs. Geologists have used the idea that a reservoir consists of the geometrically complex connected and disconnected sandbodies to be able to qualitatively describe it by means of outcrop studies and subsurface and for the reliable estimation of accessible hydrocarbons in place and the expected recoverable oil [Haldorsen et al., 1988]. Hence, sandbodies within a reservoir are assumed to be flow units connected in a complicated way through the sedimentary processes along with and poorer siltstones, mudstones and shales (of lower permeability) [Bridge and Leeder, 1979]. An example is a meandering river which deposits sand layers over the time in its flowing bed results in formation of a system of embedded sandbodies in an impermeable background (Fig. 3). Moreover, outcrop studies often show a network of connected fractures (Fig. 3). In this study, we use overlapping sandbodies and fracture network models within the framework

of continuum percolation to determine the reservoir connectivity and conductivity.

Fig. 3. Illustration of (left) Overlapping sandbodies deposited by a meandering river over

the ages (Nurafza et al 2006). (right) 2-D fracture model (Odling et al 1999)

Table 4 [Stauffer and Aharony, 1992; Sahimi, 1995].

(Chang & Yortsos 1990).

**3. Two basic models** 

�� �� ���� We assume that conductive bodies represented by squares of side size � (or circles) in two dimensional square region of side size � (or by cubes/spheres in 3D) are distributed randomly and independently in an impermeable background. The effective size of the system is defined by the dimensionless length � = � � **.** The net to gross ratio �, equivalent to the occupancy probability*,* is defined as the total area of sands divided by the total area of the region. The connectivity or connected sand fraction �(�� �)that is a function of both � and � is defined as the probability of a point on the sandbodies area, belonging to the percolating cluster. As each sandbody is put down through the model build up, it is given a color/cluster number showing a cluster that the body belongs to it. Moreover, we need an efficient computational algorithm to check the sandbody overlapping and to get the statistics of �� �� � for each realization.

#### **3.2 Fracture network model**

Many hydrocarbon reservoirs are naturally fractured. By fractures we mean any discontinuity within a rock mass that developed as a response to stress. Conductive fractures form a network of interconnected fractures with a distribution for the fracture orientation, length or aperture. The nature of fluid flow in fractured reservoirs with low matrix permeability depends strongly on the spatial distribution of the fractures. We represent fractures with line segments embedded in a 2D impermeable background (or squares/rectangles in 3D) and assume that they are the only pathways for the flow. Consider line segments representing the fractures in 2D are distributed randomly within a squared region of side size �. Fractures have the same length *l* with orientations that are distributed uniformly on the interval −90° � � � 90° . Equivalent term to the occupancy probability *p* is defined based on the probability density function that a point is in the effective area of a fracture,

$$p = 1 - e^{\frac{-N(a\_{\rm e\chi})}{4L^2}} \tag{5}$$

Where N is the number of fractures in the region and the term 〈���〉 is the average excluded area over the distribution of the fracture orientations (e.g. as shown in Masihi et al 2007 for random fracture orientation this is 〈���〉 <sup>=</sup> ��� � ). Moreover, an equivalent to the percolation probability *P* (or connected fraction) is �=����� where �� and �� are respectively the number of fractures in the spanning cluster and in the whole fracture network. Moreover, an appropriate computational algorithm is necessary to check the fracture intersections and get the statistics of �� �� � for each realizations.

Two realizations of overlapping sandbody model and fracture network model are shown in Fig. 4.

#### **4. Finite size effects**

In reality, each underground reservoir has a finite size. In such as case, instead of a single infinite cluster as occupancy increases the clusters get larger in size until one (or more) cluster spans from one side of the system to the other. For simplicity we describe this on lattice, although the argument is general and the equations describe here are applicable to other continuum models such two basic models described in the previous section. Fig. 5 shows that the connectivity from left to right may occur at very low occupancy (*p*=0.2), which is much less than the infinite percolation threshold for a 2d square lattice, but not get it at a much higher occupancy probability e.g.� = ���. As the size increases, the probability of getting such rare configurations decreases until we return to the unique single threshold value for infinite systems. In practice, for a finite system an apparent threshold ������ which depends on the system size can be used. A survey of literature [King, 1990; Stauffer and Aharony, 1992; Berkowitz, 1995; Adler and Thovert 1999; Harter, 2005] suggests that the threshold can be defined as:


Fig. 5. Connection from left to right can occur well (a) below infinite threshold or (b) above infinite threshold. (c, d) Different realizations at the same occupancy give different connectivity, *P*

We should notice that all these definitions have the same behaviour as the system size goes to infinity. The apparent threshold has the scaling [Stauffer and Aharony, 1992],

Where N is the number of fractures in the region and the term 〈���〉 is the average excluded area over the distribution of the fracture orientations (e.g. as shown in Masihi et al 2007 for

probability *P* (or connected fraction) is �=����� where �� and �� are respectively the number of fractures in the spanning cluster and in the whole fracture network. Moreover, an appropriate computational algorithm is necessary to check the fracture intersections and get

Two realizations of overlapping sandbody model and fracture network model are shown in

In reality, each underground reservoir has a finite size. In such as case, instead of a single infinite cluster as occupancy increases the clusters get larger in size until one (or more) cluster spans from one side of the system to the other. For simplicity we describe this on lattice, although the argument is general and the equations describe here are applicable to other continuum models such two basic models described in the previous section. Fig. 5 shows that the connectivity from left to right may occur at very low occupancy (*p*=0.2), which is much less than the infinite percolation threshold for a 2d square lattice, but not get it at a much higher occupancy probability e.g.� = ���. As the size increases, the probability of getting such rare configurations decreases until we return to the unique single threshold value for infinite systems. In practice, for a finite system an apparent threshold ������ which depends on the system size can be used. A survey of literature [King, 1990; Stauffer and Aharony, 1992; Berkowitz, 1995; Adler and Thovert 1999; Harter, 2005] suggests that the

� ). Moreover, an equivalent to the percolation

random fracture orientation this is 〈���〉 <sup>=</sup> ���

the statistics of �� �� � for each realizations.

Fig. 4.

**4. Finite size effects** 

threshold can be defined as:

various domain sizes.

connectivity, *P*

the point at which the connectivity is first non-zero.

the intercept of the tangent to the connectivity curve taken from its inflexion.

the point of intersection of all the curves of the probability of percolation obtained at

Fig. 5. Connection from left to right can occur well (a) below infinite threshold or (b) above

We should notice that all these definitions have the same behaviour as the system size goes

infinite threshold. (c, d) Different realizations at the same occupancy give different

to infinity. The apparent threshold has the scaling [Stauffer and Aharony, 1992],

the occupancy probability at which half of the realizations percolate.

$$
\mathfrak{p}\_c(L) - \mathfrak{p}\_c^{\infty} \propto L^{-1/v} \tag{6}
$$

From which the infinite percolation threshold �� �can be estimated. Let describe the finite size scaling for connectivity even though it is applicable to conductivity as well. For finite size systems, we define the connectivity �(�� �) as the fraction of occupied sites belonging to the spanning clusters. Again figure below shows that different realizations at the same occupancy give different values for the connectivity��(�� �).

If we plot the connected fractions��(�� �) as a function of occupancy probability "*p"* over a large number of realizations for a particular system size we get a scatter of points showing no longer a sharp transition in the connectivity behaviour near the threshold so the main effect of finite boundaries is to *smear out* the percolation transition. Fig. 6 shows a considerable spread in connectivity results for basic sandbody model with percolation threshold of �� � = 0.668 due to the finite size effects.

Fig. 6. The typical scatter plots of connectivity, *P* versus net-to-gross, *p* at system size L= 16 and 64 showing the smearing out of the percolation transition in the basic two dimensional sandbody model

Instead of working with the whole scatter we can determine the average value of connectivity �(�� �) and standard deviation of connectivity ∆(�� �) obtained over all realizations at the same occupancy probability, p. As the system size increases the amount of smearing of the percolation transition decreases and the standard deviation of connectivity curves becomes narrower and its peak becomes lower (Fig. 7).

Finite size scaling law, based on the principle that all the properties depend only on the ratio of the system size *L* to the correlation length �, is a way of relating all of the connectivity curves (i.e. *P* and ∆) obtained for different system sizes to a single curve through the scaling equations (Stauffer and Aharony, 1992):

$$P(p,L) = L^{-\beta/v} \mathfrak{S}[(p - p\_c^{\infty})L^{1/v}] \tag{7}$$

$$\Delta(p,L) = L^{-\beta/v} \Re \left[ (p - p\_c^{co}) L^{1/v} \right] \tag{8}$$

Where functions ℑ and ℜ are respectively the master function for the mean and standard deviation of connectivity results. Fig. 8 shows the data collapse presented in Figure 7 by using these scaling transformations.

Fig. 7. The average value of connectivity ���� �� and standard deviation of connectivity ���� �� results obtained over all realizations at the same occupancy probability, p at six system sizes for squared overlapping sandbody model

Fig. 8. Illustration of the finite size scaling laws for squared sandbodies in two dimensions for (a) the average connectivity and (b) standard deviation in connectivity results (Sadeghnejad et al. 2011)

These provide the master curves for the mean connected sand fraction and the standard deviation of connectivity results (ℑ and ℜ) for the case of overlapping sandbodies model from which one can predict the mean connectivity and its associated uncertainty for any system size very quickly, without performing any explicit realizations.

**Example**: Consider a typical reservoir of 50 m thick and inter-well spacing of 1 km within which sandbodies of 5 m thick and 100 m long are distributed randomely. Assuming the net-to-gross of 0.75 for the reservoir, estimate the connected sand fraction between the two wells.

(a) (b) Fig. 8. Illustration of the finite size scaling laws for squared sandbodies in two dimensions

These provide the master curves for the mean connected sand fraction and the standard deviation of connectivity results (ℑ and ℜ) for the case of overlapping sandbodies model from which one can predict the mean connectivity and its associated uncertainty for any

**Example**: Consider a typical reservoir of 50 m thick and inter-well spacing of 1 km within which sandbodies of 5 m thick and 100 m long are distributed randomely. Assuming the net-to-gross of 0.75 for the reservoir, estimate the connected sand fraction between the two

for (a) the average connectivity and (b) standard deviation in connectivity results

system size very quickly, without performing any explicit realizations.

Fig. 7. The average value of connectivity ���� �� and standard deviation of connectivity ���� �� results obtained over all realizations at the same occupancy probability, p at six

system sizes for squared overlapping sandbody model

(Sadeghnejad et al. 2011)

wells.

**Solution**: Consider a two dimensional cross section reservoir model. Then, �� <sup>=</sup> �� � = 10, ��� <sup>=</sup> ���� ��� = 10 . Moreover, we use usual values of �� � = 0.668, � � = 0.75, � � = 0.015. Then by using master curves in Figure 8 with (���� �)���� = 0.46 we get ����� = 1.05and ∆���� = 0.26 from which � = 0.82 and ∆= 0.2. Hence, the predicted connected fraction of sands between two wells is 0.82±0.20

A similar analysis can be done for the effective permeability results by considering single phase flow between two wells under a fixed pressure drop which gives the pressure equation ∇. (�∇�) = 0 where *K* is the permeability and *P* is the pressure and ∇ is the derivative operator. Having solved this on the system, one can then calculate the total flow rate and equate fluxes with those from an equivalent homogenous system to determine the effective permeability *K*. As before, once we collect the statistics of the effective permeability�(�, �) and net to gross, *p* for a given system size, we can plot the effective permeability against the net to gross that gives a scatter of points (see Fig. 9).

Fig. 9. The typical scatter plots of effective permeability � versus net-to-gross � at system size of � =16 and 64 showing the smearing out of the percolation transition in the basic two dimensional sandbody model

Again, instead of working with the whole scatter one can determine the mean effective permeability, �(�, �), and its associated uncertainty in the effective permeability results, ∆�(�, �) which has the scaling: [Stauffer and Aharony, 1992],

$$K(p,L) = L^{-\mu/\nu} \kappa[(p - p\_c^{\infty})L^{1/\nu}] \tag{9}$$

$$
\Delta\_K(p, L) = L^{-\mu/\nu} G\left[ (p - p\_c^{\infty}) L^{1/\nu} \right] \tag{10}
$$

Where functions κ and G are two master functions respectively for the mean reservoir effective permeability and its standard deviation. Fig. 10 shows how these scaling laws collapse the effective permeability results of various system size presented in Figure 9 on top of each other.

Using the master curves for the connectivity (i.e. the curve ℑ from Figure 8) and effective permeability (i.e. the curve � from Figure 10) one can replot the connectivity and effective permeability against wellspacing at various net to gross values (Fig. 11). The connected sand fraction shows the fraction of the original oil in place that is in contact with two wells.

Fig. 10. Illustration of (left) effective permeability results �(�� �)of squared sandbodies in two dimensions at various system sizes (right) the data collapse for conductivity results using finite size scaling

Fig. 11. Illustration of the connectivity results�(�� �)and effective permeability results �(�� �)of squared overlapping sandbodies against well spacing L at various net-to-gross *p* values

#### **5. Anisotropic reservoirs**

In isotropic reservoirs, the horizontal connectivity is the same as the vertical connectivity on average if not for individual realizations. However, for many realistic systems, the objects or their orientation is rarely isotropic. For example in fractured rocks, fracture sets with particular orientations are typically formed as a result of tectonic history [Bear *et al*., 1993; Sahimi, 1995; Adler and Thovert, 1999; Zhang and Sanderson, 2002] or in overlapping sandbody model the simple assumption of squared sands may not be true and rectangles can be used. This leads to the creation of an *easy* direction for connected paths which is in

Fig. 10. Illustration of (left) effective permeability results �(�� �)of squared sandbodies in two dimensions at various system sizes (right) the data collapse for conductivity results

In isotropic reservoirs, the horizontal connectivity is the same as the vertical connectivity on average if not for individual realizations. However, for many realistic systems, the objects or their orientation is rarely isotropic. For example in fractured rocks, fracture sets with particular orientations are typically formed as a result of tectonic history [Bear *et al*., 1993; Sahimi, 1995; Adler and Thovert, 1999; Zhang and Sanderson, 2002] or in overlapping sandbody model the simple assumption of squared sands may not be true and rectangles can be used. This leads to the creation of an *easy* direction for connected paths which is in

Fig. 11. Illustration of the connectivity results�(�� �)and effective permeability results �(�� �)of squared overlapping sandbodies against well spacing L at various net-to-gross *p*

using finite size scaling

values

**5. Anisotropic reservoirs** 

the short direction and a *difficult* direction which is along the long axis. Monte-Carlo, conformal field and renormalization group theory and duality arguments can be used to incorporate the anisotropic effects in the finite size scaling laws. We define an *apparent* percolation threshold in each direction to be the point where 50% of realizations connect in that direction. Again for simplicity, we describe this on lattice although the results are applicable to other continuum systems such as overlapping sandbody system or fracture system (Masihi et al 2006, 2007; Sadeghnejad et al 2010). Consider an anisotropic system in 2d space which now has different effective size *Lx* and *Ly* along the x and y directions and define an aspect ratio *ω = Lx/Ly* to represent the anisotropy. Obviously if the effective sizes in two directions are equal (i.e. *Lx = Ly*) then *ω =1* which means that the system is isotropic (see Fig. 12 for a typical realization).

Fig. 12. Illustration of three connected clusters in the Y direction at *ω =10* and *p=0.55* on a lattice which is globally anisotropic (*Lx =200, Ly=20*)

Fig. 13. Illustration of threshold in the horizontal direction for overlapping sandbody model which is globally isotropic but the objects are rectangles with the number of objects N=280, the net to gross *p*=0.659 and the horizontal connected sand fraction of *Ph*=0.449

We investigate finite size scaling at a fixed aspect ratio�� and with the length scale *Lx*. We expect that both of the apparent *x* and *y* percolation thresholds to scale with the length scale as,

$$
\mathfrak{p}\_c^l - \mathfrak{p}\_c^{\infty} = \Lambda\_l(\omega) L\_\chi^{-1/\nu} \tag{11}
$$

Where ��� � is the apparent threshold and Λ� is simply the constant of proportionality which is dependent on the aspect ratio ω with *i* labels the coordinate direction. The following simple scaling function has been observed for the proportionality constant Λ�:

$$\Lambda(\omega) = \mathbf{c} \{ \omega^{1/\nu} - \mathbf{1} \} \tag{12}$$

Where the coefficient *c* has been estimated numerically about 0.92, 0.45 and 0.58 respectively for the lattice model, overlapping sandbody model and fracture model (King 1990, Masihi et al 2006, 2007, Sadeghnejad et al 2011). It was found that using infinite percolation threshold within the usual finite size scaling law then the anisotropic connectivity in the horizontal and vertical directions (i.e. �� and ��) collapse onto the single universal isotropic connectivity curve (ℑ) as shown in Fig. 14.

Fig. 14. Illustration of finite size scaling for the average connected fraction (�) of overlapping squares in horizontal and vertical directions by using infinite threshold (left) and apparent threshold (right) in the finite size scaling law.

This led to a way to modify the finite size scaling in the presence of anisotropy. Such scaling, for example for connectivity results becomes (Masihi et al 2006, 2007):

$$P(p, L\_{\chi}, \omega) = L\_{\chi}^{-\beta/\nu} \mathfrak{S} \Big[ (p - \mathfrak{p}\_c) L\_{\chi}^{-1/\nu} \Big] \tag{13a}$$

$$\Delta(p, L\_{\chi}, \omega) = \omega^{1/2} L\_{\chi}^{-\beta/v} \Re \left[ (p - \mathfrak{p}\_c) L\_{\chi}^{-1/v} \right] \tag{13b}$$

These results enable us to use the same isotropic universal curves given in Figure 8 for predicting the connectivity of anisotropic lattices which is a good enough approximation for engineering purposes.

**Example**: Plot the variation of the connected sand fraction against well spacing at net-togross values of *p=*0.65, 0.7, 0.75 and 0.80 for a cross section reservoir model of 10m thickness which contains sandbodies of 50 m long and 1 m thick.

**Solution**: with a well spacing of 500 m then the effective system size is �� = 10, �� = 10 so the system is isotropic. To calculate the connected sand fraction at a particular net-to-gross (*p*) we need to invert the finite size scaling with *L*=10 and the master curve in Fig.8 At a well spacing of 1000 m the effective system size is �� = 20, �� = 10 so the system is anisotropic with aspect ratio of ω=2. The horizontal connectivity of this system is the same as the vertical connectivity of a system with � = 10 and the aspect ratio of ω = 0.5. With this we can determine the shift in apparent threshold by using Eq. 12 which is Λ(ω) = −0.18 and so the apparent threshold from Eq. 11 becomes��� = 0.668 + 0.18 × 10��.�� = 0.7.Again the

Where the coefficient *c* has been estimated numerically about 0.92, 0.45 and 0.58 respectively for the lattice model, overlapping sandbody model and fracture model (King 1990, Masihi et al 2006, 2007, Sadeghnejad et al 2011). It was found that using infinite percolation threshold within the usual finite size scaling law then the anisotropic connectivity in the horizontal and vertical directions (i.e. �� and ��) collapse onto the single universal isotropic

 Fig. 14. Illustration of finite size scaling for the average connected fraction (�) of overlapping squares in horizontal and vertical directions by using infinite threshold (left) and apparent

This led to a way to modify the finite size scaling in the presence of anisotropy. Such scaling,

These results enable us to use the same isotropic universal curves given in Figure 8 for predicting the connectivity of anisotropic lattices which is a good enough approximation for

**Example**: Plot the variation of the connected sand fraction against well spacing at net-togross values of *p=*0.65, 0.7, 0.75 and 0.80 for a cross section reservoir model of 10m thickness

**Solution**: with a well spacing of 500 m then the effective system size is �� = 10, �� = 10 so the system is isotropic. To calculate the connected sand fraction at a particular net-to-gross (*p*) we need to invert the finite size scaling with *L*=10 and the master curve in Fig.8 At a well spacing of 1000 m the effective system size is �� = 20, �� = 10 so the system is anisotropic with aspect ratio of ω=2. The horizontal connectivity of this system is the same as the vertical connectivity of a system with � = 10 and the aspect ratio of ω = 0.5. With this we can determine the shift in apparent threshold by using Eq. 12 which is Λ(ω) = −0.18 and so the apparent threshold from Eq. 11 becomes��� = 0.668 + 0.18 × 10��.�� = 0.7.Again the

����ℑ�(�−���)��

������(�−���)��

����� (13a)

����� (13b)

for example for connectivity results becomes (Masihi et al 2006, 2007):

�(�, ��, ω) = ��

∆(�, ��, ω) = ω�����

connectivity curve (ℑ) as shown in Fig. 14.

threshold (right) in the finite size scaling law.

which contains sandbodies of 50 m long and 1 m thick.

engineering purposes.

connected sand fraction at a particular net-to-gross (*p*) can be calculated by inverting the finite size scaling with *L*=10 and by using the numerical values from the connectivity master curve. Now by continuing this process we can determine the connectivity as a function of well spacing for various net-to-gross ratios (see Fig. 11 for typical behavior).

#### **6. Size variation and orientation distribution**

In reality, the sandbodies in the overlapping sandbody model or fractures in the fracture network model may not have the same size or they do not have a fixed orientation. The distributions for the size and orientation are related to the sedimentological environment which deposited the reservoir. The finite size scaling discussed so far assumes that the sandbodies all have the same lengths and the global behaviour is controlled by two lengths, the system size, *L*, and the correlation length, *ξ.* If there is distribution of sizes then this is changed. If the distribution of sizes is a continuous (e.g. uniform or a Gaussian distribution), then it was suggested to replace this variable size system with another system in which the bodies are all with the same effective size. This means that the connectivity of, for example, sandbodies of variable size is identical to the connectivity of sandbodies of the same size. Consider the sandbodies that are represented by squares of variable side size of �. Then it was shown that the effective size ���������� which gives the same percolation behavior can be based on the second moment of size distribution as[Balberg et al 1984, Masihi et al 2008]:

$$l\_{effective} = \sqrt{\langle a^2 \rangle} \tag{14}$$

 Where 〈 〉 denotes the average over all fracture lengths. This hypothesis has been verified on the 2d overlapping sandbody model [King 1990], two and three dimensional fracture model [Masihi et al 2008] and 3d overlapping sandbody model [Sadeghnezad et al 2011]. Fig 15 shows the connectivity and conductivity results for a system with size variation as compared with the master curves when effective size is used.

Fig. 15. (left) connected sand fractions and (right) effective permeability for variable sizes squares with uniform length distribution in the range of [0.5-1.5] in a system of size L=24

$$p\_{\mathcal{C}} = 1 - e^{-\rho\_{\mathcal{C}}} \tag{15}$$

$$\rho\_c = \frac{4.41}{2\left(1 + \frac{\sin^2 \theta\_o}{\theta\_0^2}\right) + \frac{1 + \omega^2}{2\omega \theta\_0^2} (2\theta\_o - \sin 2\theta\_o)}\tag{16}$$

$$
\langle \ell\_x \rangle = \frac{1}{n} \Sigma \{a \cos \theta + b \sin |\theta|\}, \\
\langle \ell\_y \rangle = \frac{1}{n} \Sigma (b \cos \theta + a \sin |\theta|) \tag{17}
$$

al., 1999; Andrade et al. 2000; King et al., 2002) as well as the post breakthrough behavior. The knowledge of these is important, for example, in water flooding to improve recovery in hydrocarbon reservoirs. The conventional approach to predict the breakthrough time or post breakthrough behavior is to build a detailed geological model, upscale it, and then perform flow simulations. To estimate the uncertainty in such a system, this has to be repeated for various realizations of the geological model. The problem with this approach is that it is computationally very expensive and time consuming. Although there has been some progress to reduce the calculation time [Donato and Blunt, 2004]; there is a great incentive to produce much simpler physically-based techniques to predict these reservoir performance parameters very quickly, particularly for engineering purposes. In this section, we use percolation-based reservoir structures [Stauffer and Aharony 1994] on which we can analyze the breakthrough time and post breakthrough behavior between a injector and a producer [Andrade et al 2000].

Dokholyan et al. [1999] did the pre-elementary study of the breakthrough time. They presented a scaling ansatz for the distribution function of the shortest paths connecting any two points on a percolating cluster. Traveling time and traveling length for tracer dispersion in two-dimensional bond percolation systems have been studied by Lee et al. [1999]. King et al. [1999] used percolation theory to predict the distribution of the shortest path between well pairs and presentd a scaling hypothesis for this distribution which has been confirmed by the numerical simulation. Andrade et al. [2000] concentrated on the flow of fluid between two sites on the percolation cluster. They modelled the medium by using bond percolation on a lattice, while the flow front was modeled by tracer particles driven by a pressure difference between two fixed sites representing the injection and productions wells. They investigated the distribution function of the shortest path connecting these two wells, and proposed a scaling which was confirmed by extensive simulations. Moreover, Andrade et al. [2001] investigated the dynamics of viscous penetration in two-dimensional percolation networks at criticality for the case in which the ratio between the viscosities of displaced and injected fluids is very large. They reported extensive numerical simulations showing that the scaling exponents for the breakthrough time distribution was the same as the previously reported values computed for the case of unit viscosity ratio. Araujo et al. [2002] analysed the distributions of traveling length and minimal traveling time through two-dimensional percolation porous media characterized by relatively long-range spatial correlations. They found that the probability distribution functions follow the same scaling ansatz originally proposed for the uncorrelated case, but with quite different scaling exponents. Using Monte Carlo simulations, comparing the results of the proposed scaling laws from percolation theory with the time for a fluid injected into an oil field to breakthrough into a production well was studied by King et al. [2002]. Lopez et al. [2003] numerically simulated the traveling time of a tracer in convective flow between two wells in a percolation porous media. They analyzed the traveling time probability density function for two values of the fraction of connecting bonds, namely the homogeneous case and the inhomogeneous critical threshold case. Soares et al. [2004] investigated the distribution of shortest paths in percolation systems at the percolation threshold from one injector well to multiple producer wells. They calculated the probability distribution of the optimal path length. Li et al. [2009] applied percolation method to estimate inter-well connectivity of thin intervals for noncommunicating stratigraphic intervals in an oil field.

al., 1999; Andrade et al. 2000; King et al., 2002) as well as the post breakthrough behavior. The knowledge of these is important, for example, in water flooding to improve recovery in hydrocarbon reservoirs. The conventional approach to predict the breakthrough time or post breakthrough behavior is to build a detailed geological model, upscale it, and then perform flow simulations. To estimate the uncertainty in such a system, this has to be repeated for various realizations of the geological model. The problem with this approach is that it is computationally very expensive and time consuming. Although there has been some progress to reduce the calculation time [Donato and Blunt, 2004]; there is a great incentive to produce much simpler physically-based techniques to predict these reservoir performance parameters very quickly, particularly for engineering purposes. In this section, we use percolation-based reservoir structures [Stauffer and Aharony 1994] on which we can analyze the breakthrough time and post breakthrough behavior between a injector and a

Dokholyan et al. [1999] did the pre-elementary study of the breakthrough time. They presented a scaling ansatz for the distribution function of the shortest paths connecting any two points on a percolating cluster. Traveling time and traveling length for tracer dispersion in two-dimensional bond percolation systems have been studied by Lee et al. [1999]. King et al. [1999] used percolation theory to predict the distribution of the shortest path between well pairs and presentd a scaling hypothesis for this distribution which has been confirmed by the numerical simulation. Andrade et al. [2000] concentrated on the flow of fluid between two sites on the percolation cluster. They modelled the medium by using bond percolation on a lattice, while the flow front was modeled by tracer particles driven by a pressure difference between two fixed sites representing the injection and productions wells. They investigated the distribution function of the shortest path connecting these two wells, and proposed a scaling which was confirmed by extensive simulations. Moreover, Andrade et al. [2001] investigated the dynamics of viscous penetration in two-dimensional percolation networks at criticality for the case in which the ratio between the viscosities of displaced and injected fluids is very large. They reported extensive numerical simulations showing that the scaling exponents for the breakthrough time distribution was the same as the previously reported values computed for the case of unit viscosity ratio. Araujo et al. [2002] analysed the distributions of traveling length and minimal traveling time through two-dimensional percolation porous media characterized by relatively long-range spatial correlations. They found that the probability distribution functions follow the same scaling ansatz originally proposed for the uncorrelated case, but with quite different scaling exponents. Using Monte Carlo simulations, comparing the results of the proposed scaling laws from percolation theory with the time for a fluid injected into an oil field to breakthrough into a production well was studied by King et al. [2002]. Lopez et al. [2003] numerically simulated the traveling time of a tracer in convective flow between two wells in a percolation porous media. They analyzed the traveling time probability density function for two values of the fraction of connecting bonds, namely the homogeneous case and the inhomogeneous critical threshold case. Soares et al. [2004] investigated the distribution of shortest paths in percolation systems at the percolation threshold from one injector well to multiple producer wells. They calculated the probability distribution of the optimal path length. Li et al. [2009] applied percolation method to estimate inter-well connectivity of thin intervals for non-

producer [Andrade et al 2000].

communicating stratigraphic intervals in an oil field.

In this part, we characterize the breakthrough time and post breakthrough behaviour between an injector and a producer in a real case by using the Burgan reservoir dataset of Norouz offshore oil field in the south of Iran and applying the propsoed scaling law against the detailed reservoir simulation results. Furthermore, we develop a scaling law for the prediction of post breakthrough behavior based on the probability distribution of breakthrough time. For this purpose, we focus on the time to reach to the water cut of 50% in the production well and introduce a new probability distribution function for this.

It has been shown that the distribution of breakthrough time ��� between two wells with well spacing � conditioned to the formation size, *L*, and net-to-gross ratio, *p*, has the scaling [Andrade et al 2000]:

$$P(t\_{br}|r,L,p) \propto \frac{1}{r^{d\_t}} \left(\frac{t\_{br}}{r^{d\_t}}\right)^{-g\_t} f\_1\left(\frac{t\_{br}}{r^{d\_t}}\right) f\_2\left(\frac{t\_{br}}{L^{d\_t}}\right) f\_3\left(\frac{t\_{br}}{\xi^{d\_t}}\right) \tag{18}$$

Where functions ��(�) � ������ , ��(�) � ������ and ��(�) � ������ . Moreover, ���is the breakthrough time for a given realization, *r* is the Euclidean distance between the injector and a producer, and � is the typical size of reservoir sandbodies. The numerical values of the other exponents were determined through rigorous simulation studies (Table 5).


Table 5. Numerical values of the exponent in Eq. 18[Andrade et al 2000]

Eq. 18 assumed that the injecting and producing fluids are incompressible, the displacing fluid has the same viscosity and density as displaced fluid (i.e. like passive tracer transport). Pressure field is determined by the solution of the single phase flow equation (i.e. Laplace equation, �� (���) � �. Fixed pressure boundary conditions were considered at the two wells. Breakthrough time corresponds to the time when the first tracer reaches the production well for a given reservoir realization. The probability function given by Eq. 1 can be encoded in a spreadsheet from which, using some primary data, the probability distribution of the breakthrough time is determined very quickly. The interested reader is referred to work of Andrade et al 2000 for the idea used in developing Eq. 18

#### **7.1 Application to real field**

To validate the approach, we have used the Burgan formation dataset of Norouz offshore oil field in the south of Iran. Core and palynological data have indicated that the Burgan consists of a series of incision-fill sequences occurring in an estuarine/coastal plain/deltaic environment. Consequently, the Burgan formation consists of a thick stack of excellent quality sands incising into each other with a few remaining shalier sediments locally separating these sequences. The excellent formation sands of the Burgan consist mainly of valley-fill deposition sediments, where freshwater fluvial, coarse-grained sediment (high permeable zones) accumulated in the upper and middle reaches of the incised valley system, while tidally influenced sediments accumulated in the marine-influenced middle and lower reaches. As a result of this aggradations, less fluvial-derived sediment reached the lower reaches of the valley system, which resulted in the clean, well-sorted shoreline material being re-worked and deposited in the tidal channels in this area [Huerlimann 2004]. Because of this process, the valley-fill deposition sediments can be considered as high permeable flow units in a background that is essentially impermeable which makes the Burgan reservoir ready to be modeled with the percolation approach.

To have different realizations, various points for the well locations in the entire Burgan reservoir were randomly selected. Nearly 300 flow simulations were run on these wells at the full field scale of the Burgan reservoir formation. To include well spacing effect, two distinct well spacing has been considered as 500m, and 1000m. Then we plot the histogram of the breakthrough times by using the results of 300 simulation runs. Moreover, we encoded the Eq. 18 in an excel spreadsheet for the Burgan reservoir data to produce the breakthrough time probability distribution. The detail of this is as follows: We first make the reservoir sizes dimensionless with the dimension of the sandbody in the appropriate direction (��� ��� ��) and apply the scaling for the ��� with the minimum dimensionless length involved ���� = ������� ��� ��). We note that the scaling was for the passive tracer transport; however, the actual displacement is not the same as passive tracer transport. Hence, the Buckley Leverett displacement idea along fastest streamline may be assumed. Moreover, it should be noted that in encoding Eq. 18 with the scaling ����� we actually mean � �� � � � �� � �� where �� is a typical length in the system (e.g. sandbody size or ��) and �� is the typical time (e.g. the time needed to transit one sandbody). Consider the transit time � (in seconds) for a fluid of viscosity � (cp) between two wells of radius ��(in cm) with pressure drop �� (in atm) separated by a distance � (in cm) in a homogeneous reservoir with a permeability � (Darcy) is given by � = ���� ���� �� � �� . Then, with linear assumption for the pressure drop across the well spacing, the pressure drop across the sandbody becomes �� � � �� �. Taking into account the effect of the sandbodies overlapping through the simulation analysis then the typical transit time for a sandbody becomes, �� <sup>=</sup> ����� � ������ � �� � �� �� �� . These procedures were repeated for two cases with well spacings 500 and 1000 m. The comparison of the result of the scaling law of the breakthrough time (curve) with the results obtained from the conventional numerical simulations (bar chart) for two well spacing is shown in Fig. 19.

To validate the approach, we have used the Burgan formation dataset of Norouz offshore oil field in the south of Iran. Core and palynological data have indicated that the Burgan consists of a series of incision-fill sequences occurring in an estuarine/coastal plain/deltaic environment. Consequently, the Burgan formation consists of a thick stack of excellent quality sands incising into each other with a few remaining shalier sediments locally separating these sequences. The excellent formation sands of the Burgan consist mainly of valley-fill deposition sediments, where freshwater fluvial, coarse-grained sediment (high permeable zones) accumulated in the upper and middle reaches of the incised valley system, while tidally influenced sediments accumulated in the marine-influenced middle and lower reaches. As a result of this aggradations, less fluvial-derived sediment reached the lower reaches of the valley system, which resulted in the clean, well-sorted shoreline material being re-worked and deposited in the tidal channels in this area [Huerlimann 2004]. Because of this process, the valley-fill deposition sediments can be considered as high permeable flow units in a background that is essentially impermeable which makes the

To have different realizations, various points for the well locations in the entire Burgan reservoir were randomly selected. Nearly 300 flow simulations were run on these wells at the full field scale of the Burgan reservoir formation. To include well spacing effect, two distinct well spacing has been considered as 500m, and 1000m. Then we plot the histogram of the breakthrough times by using the results of 300 simulation runs. Moreover, we encoded the Eq. 18 in an excel spreadsheet for the Burgan reservoir data to produce the breakthrough time probability distribution. The detail of this is as follows: We first make the reservoir sizes dimensionless with the dimension of the sandbody in the appropriate direction (��� ��� ��) and apply the scaling for the ��� with the minimum dimensionless length involved ���� = ������� ��� ��). We note that the scaling was for the passive tracer transport; however, the actual displacement is not the same as passive tracer transport. Hence, the Buckley Leverett displacement idea along fastest streamline may be assumed. Moreover, it should be noted that in encoding Eq. 18 with the scaling ����� we actually

where �� is a typical length in the system (e.g. sandbody size or ��) and �� is

. Then, with linear assumption for the

� ������ � �� � �� �� ��

. These

the typical time (e.g. the time needed to transit one sandbody). Consider the transit time � (in seconds) for a fluid of viscosity � (cp) between two wells of radius ��(in cm) with pressure drop �� (in atm) separated by a distance � (in cm) in a homogeneous reservoir

> ���� �� � ��

�. Taking into account the effect of the sandbodies overlapping through the simulation

pressure drop across the well spacing, the pressure drop across the sandbody becomes

procedures were repeated for two cases with well spacings 500 and 1000 m. The comparison of the result of the scaling law of the breakthrough time (curve) with the results obtained from the conventional numerical simulations (bar chart) for two well spacing is shown in

analysis then the typical transit time for a sandbody becomes, �� <sup>=</sup> �����

Burgan reservoir ready to be modeled with the percolation approach.

with a permeability � (Darcy) is given by � = ����

**7.1 Application to real field** 

mean � �� � � � �� � ��

�� � � ��

Fig. 19.

Fig. 19. Comparison of prediction of the breakthrough time using the scaling laws of Eq. 18 (curve) with results from the numerical simulations (bar chart) for well spacing of a) 500m, b) 1000 m.

As it can be seen from Fig. 19, there is a reasonable agreement between the prediction from the scaling law, given by Eq. 18, and the direct numerical simulation results.

#### **7.2 Post breakthrough behavior**

The main question in evaluating the post breakthrough behavior is that if the probability distribution of the breakthrough time is available how does this reduce the uncertainty in the post breakthrough behavior. Specifically, we want to check if there is a correlation between the breakthrough time results and the time taken for the oil production to fall by, for example 50%, so called ݐଵȀଶ, (or water cut to increase to 50%). To get these statistics, we continue the conventional flow simulations to reach to this water cut value. Then we plot ݐଵȀଶ versus ݐ on a log-log scale. Fig. 20 shows the results these cross plots at two well spacing values.

Fig. 20. Illustration of time to reach to the water cut of 50% in production wells versus breakthrough time for different well spacing as a)500m, b) 1000 m.

As it can be seen, there is a positive correlation between the breakthrough time and the time to reach to the water cut of 50% in the production well. The degree of correlation observed in Fig. 20 depends mainly on the heterogeneity of the reservoir. It should be noted that the understudy reservoir is a heterogeneous channel reservoir. By plotting the ݐଵȀଶ versus ݐ of all simulation results from different well spacing on a log-log scale, one can find a power law scaling between the two times ݐଵȀଶ and ݐ as:

$$t\_{1/2} = 9.343t\_{br}^{1.27} \tag{19}$$

Therefore, for estimation of the probability distribution of the post breakthrough time (defined by the water cut of 50%), one can use the probability distribution of the breakthrough time from Eq.18 and scaled it by using Eq. 19. Fig. 21 compares the simulation results of the post breakthrough behavior with the proposed scaling law.

Fig. 21. Comparison of prediction of the time to reach to the water cut of 50% using the scaling laws of post breakthrough behavior (curve) with results from the numerical simulations (bar chart)

#### **8. Summary and conclusion**

We have described an approach based on percolation theory to model underground reservoirs and to estimate the reservoir performance parameters very rapidly. This is a field scale approach and it uses simple data acquired from reservoir geology such as reservoir size, netto-gross ratio within the framework of continuum percolation theory to rapidly evaluate the important parameters of geometrically complex reservoirs, especially in the early reservoir development stage when data are very limited. This approach has been used in estimating the connectivity and conductivity between two or more wells in the reservoir. These can be used for reserve estimation and infill drilling projects. Moreover much details of the methodology have been presented. In particular, the effects of anisotropy in the system, size variations and angular dispersion of sandbodies on the percolation results have been addressed. We observed that the effective length based on the second moment of length distribution can be used to incorporate the effect of variable size. Moreover, the approach to incorporate the effects of angular dispersion of sandbodies within the finite size scaling has been presented.

The methodology has not been limited to estimate the reservoir static data. We have described how to estimate the breakthrough time and post breakthrough behavior between two injection and production wells in terms of percolation concepts. In particular, the proposed scaling law of the breakthrough time probability has been compared with the results obtained from the conventional numerical simulations on the Burgan formation dataset of Norouz offshore oil field in the south of Iran which showed a good agreement between the prediction from the percolation based expression and the numerical simulation results. We have also extend the approach to estimate the post breakthrough behavior (e.g. the time to reach to the water cut of 50% in the production well) given the breakthrough time is available. There was a good agreement between the prediction from the derived correlation and the numerical simulation results. Moreover, the prediction from the scaling law took a fraction of a second of CPU times (as it only needs some algebraic calculations) compared with many hours required for the conventional numerical simulations.

#### **9. References**

284 Water Resources Management and Modeling

all simulation results from different well spacing on a log-log scale, one can find a power

ݐଵȀଶ ൌ ͻǤ͵Ͷ͵ݐ

Fig. 21. Comparison of prediction of the time to reach to the water cut of 50% using the scaling laws of post breakthrough behavior (curve) with results from the numerical

We have described an approach based on percolation theory to model underground reservoirs and to estimate the reservoir performance parameters very rapidly. This is a field scale approach and it uses simple data acquired from reservoir geology such as reservoir size, netto-gross ratio within the framework of continuum percolation theory to rapidly evaluate the important parameters of geometrically complex reservoirs, especially in the early reservoir development stage when data are very limited. This approach has been used in estimating the connectivity and conductivity between two or more wells in the reservoir. These can be used for reserve estimation and infill drilling projects. Moreover much details of the methodology have been presented. In particular, the effects of anisotropy in the system, size variations and angular dispersion of sandbodies on the percolation results have been addressed. We observed that the effective length based on the second moment of length distribution can be used to incorporate the effect of variable size. Moreover, the approach to incorporate the effects of

angular dispersion of sandbodies within the finite size scaling has been presented.

The methodology has not been limited to estimate the reservoir static data. We have described how to estimate the breakthrough time and post breakthrough behavior between two injection and production wells in terms of percolation concepts. In particular, the proposed scaling law of the breakthrough time probability has been compared with the results obtained from the conventional numerical simulations on the Burgan formation dataset of Norouz offshore oil field in the south of Iran which showed a good agreement between the prediction from the percolation based expression and the numerical simulation

results of the post breakthrough behavior with the proposed scaling law.

Therefore, for estimation of the probability distribution of the post breakthrough time (defined by the water cut of 50%), one can use the probability distribution of the breakthrough time from Eq.18 and scaled it by using Eq. 19. Fig. 21 compares the simulation

ଵǤଶ

(19)

law scaling between the two times ݐଵȀଶ and ݐ as:

simulations (bar chart)

**8. Summary and conclusion** 

