Percolation Approach in Underground Reservoir Modeling

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.


Introduction
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).

www.intechopen.com
Water Resources Management and Modeling 264 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. 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.

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 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. 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 ∞ 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 andBalberg, 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]: where p ∞ 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], With as the correlation length exponent. Moreover, the mean cluster size follows a simple power law as, Where the is called the mean cluster size exponent. 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]: 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). 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 e x p o n e n t s a s i n t h e c a s e o f l a t t i c e p e r c o l a t i o n . H o w e v e r , w e n e e d t o t h e p e r c o l a t i o n 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 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 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).

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.

Overlapping sandbody model
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.

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.

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:  the point at which the connectivity is first non-zero.  the intercept of the tangent to the connectivity curve taken from its inflexion.  the occupancy probability at which half of the realizations percolate.  the point of intersection of all the curves of the probability of percolation obtained at various domain sizes. 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], 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 = .
due to the finite size effects. 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): 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. 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.
from which = . and ∆= . . 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 ∇. ∇ = 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). , , and its associated uncertainty in the effective permeability results, ∆ , which has the scaling: [Stauffer and Aharony, 1992], 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

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 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 , 2007Sadeghnejad et al 2010). Consider an anisotropic system in 2d space which now has different effective size L x and L y along the x and y directions and define an aspect ratio ω = L x /L y to represent the anisotropy. Obviously if the effective sizes in two directions are equal (i.e. L x = L y ) 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 (L x =200, L y =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 P h =0.449 We investigate finite size scaling at a fixed aspect ratio and with the length scale L x . We expect that both of the apparent x and y percolation thresholds to scale with the length scale as, 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 Λ : 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, 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. 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 , 2007: 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 = , = 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 = , = so the system is anisotropic with aspect ratio of ω= . The horizontal connectivity of this system is the same as the vertical connectivity of a system with = and the aspect ratio of ω = . . With this we can determine the shift in apparent threshold by using Eq. 12 which is Λ ω = − . and so the apparent threshold from Eq. 11 becomes = .
+ . × 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).

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: 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. Moreover, the orientational disorder of the objects (sands or fractures) may enhance the percolation (e.g. connectivity) which means that the infinite threshold has been reduced. In addition to this shift in the threshold there will be a finite size "effective" shift. Consider a system of sandbodies with uniform orientation between − , . To account for the shift in the finite threshold a conjecture about the universality of percolation thresholds [Balberg et al 1984, Balberg 1985 saying that there is a universal value of the excluded area (i.e. the area around a geometrical object in which the centre of another object must lie in order for them to overlap) of percolating shapes at the threshold was used [King 1990, Masihi et al 2008 which led to the following infinite threshold: For example, in an overlapping sandbody model with rectangles of an aspect ratio and a uniform distribution between − , the reduction of the percolation threshold as a function of aspect ratio becomes [ Fig. 16], that is in agreement with the numerical results.
To account for the finite size "effective" shift due to the orientational disorder, we use the concept of the effective length which gives the total "reach" of a body. In Fig 18 the  Where sum is over all sandbodies and n is the number of sandbodies within the rectangular system of size and . This shows that the effect of the orientational disorder is to make the sandbodies a bit larger and a bit less elongated. Then, the effective size of the system is defined by two dimensionless lengths = / ℓ and = / ℓ which implies a new aspect ratio, , as = / and so affect the finite size percolation threshold (Eq. 11) Fig. 18. Illustration of the effective size of a rectangular sandbody of angle in two directions Extensions of the methodology presented in this chapter to three dimensions can be found elsewhere (Masihi et al 2007, Sadeghnejad et al 2010. Moreover, a number of case study including using the data set of a carbonate gas condensate reservoir with tilted fault blocks in the Gulf of Gabes, offshore Tunisia , a fractured limestone bed exposed on the southern margin of the Bristol Channel Basin the north Somerset, UK [Masihi et al 2007], and the Burgen reservoir dataset of Nowrooz offshore oil field in the north of Persian Gulf, Iran [Sadeghnejad et al 2011] has been investigated which proved the applicability of percolation approach.

Breakthrough and post breakthrough behavior
Percolation concepts can be used to find the probability distribution for the breakthrough time t br between an injector and a producer (Dokholyan et al., 1998;King et 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].  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].  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.
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]: 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)  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

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 =min , , . 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, = ∆ . 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.
(a) (b) 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.

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. 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: 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)

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.