Fractal Growth Model for Estimating Breakthrough Time and Sweep Efficiency When Waterflooding Geologically Heterogeneous Rocks

We describe a fast method for estimating flow through a porous medium with a heterogeneous permeability distribution. The main application is to contaminant transport in aquifers and recovery of oil by waterflooding, where such geological heterogeneities can result in regions of bypassed contaminants or oil. The extent of this bypassing is normally assessed by a numerical flow simulation that can take many hours of computer time. Ideally the impact of uncertainty in the geological description is then evaluated by the performing of many such simulations using different realizations of the permeability distribution. Obviously, a proper Monte Carlo evaluation may be impossible when the flow simulations are so computationally intensive. Consequently, methods from statistical mechanics, such as percolation theory and random walkers (such as diffusion-limited aggregation), have been proposed; however, these methods are limited to geological heterogeneities where the correlation lengths are smaller than the system size or to continuous permeability distributions. Here we describe a growth model that can be used to estimate the breakthrough time of the water (and hence the sweep efficiency) in most types of geologically heterogeneous rocks. We show how the model gives good estimates of the breakthrough time of water at the production well in a fraction of the time needed to perform a full flow simulation.


I. INTRODUCTION
Single-phase flow and multiphase flow through porous systems play a key role in many scientific and technological areas, such as hydrogeology and petroleum engineering. In those systems geological heterogeneity can affect the flow of water, creating regions of bypassed contaminants or oil. The time of water breakthrough and the extent of these bypassed zones are normally assessed by numerical models that solve the associated partial differential equations (mass balance and the multifluid extension of Darcy's law) using a finite-volume approach, but each such simulation can take many hours of computer time. A further problem with this approach is that the exact pattern of hydraulic conductivity in the subsurface is usually uncertain and yet a very fine grid is needed to resolve all the different length scales of geological heterogeneity that may impact the flow. The impact of uncertainty is usually assessed by a Monte Carlo simulation: numerous realizations of the hydraulic conductivity pattern are created and then an expensive, multiphase simulation of flow is performed for each realization. This can be computationally impractical.
For these reasons a range of methods from statistical mechanics, such as percolation [1][2][3][4][5], invasion percolation * p.gago@imperial.ac.uk [6,7], continuum random walkers [8,9], and diffusionlimited aggregation [10], porous media. These models have been shown to be potentially very useful for estimating both the breakthrough time of the displacing fluid and some of the bulk properties of the heterogeneous system, such as effective permeability. They are especially useful for systems with a statistically homogeneous permeability distribution or when the spatial permeability distribution has a short correlation length compared with the system size. Unfortunately, these methods are less applicable to geological systems that have spatial correlation lengths that are similar to the system size.
In this work we introduce a variation of Eden's growth model [11] to estimate displacing-fluid breakthrough time in heterogeneous porous media with spatial correlation lengths similar to the system size. Eden's growth model is one type of fractal growth model and falls into the group of models classified as cellular automata. In Eden's model, growth moves outward from existing occupied sites with a probability related to the number of occupied neighbors. In our model, the probability of growth outward from an occupied site is calculated from the permeability of that occupied site. The model is different from diffusion-limited aggregation [10,12], where growth occurs by points diffusing toward the growing cluster from outside the cluster and then sticking to the surface. While diffusion-limited aggregation was used to estimate sweep efficiency due to viscous fingers [10] assuming an infinite viscosity ratio between the fluids, here we assume a mobility ratio of 1. This application focuses on flow at length scales of several meters so that capillary effects can be ignored. We demonstrate our model by considering the case of injection of water into a heterogeneous oil reservoir with the aim of rapidly estimating the dimensionless breakthrough time of the water and hence the sweep efficiency.
This work is structured as follow. Section II is divided into two parts: the first part introduces the growth model with its rules and assumptions; in the second part we review some relevant concepts of multiphase flow in porous media necessary to relate cluster growth in a grid to breakthrough time in an oil reservoir. In Sec. III we describe the geological models that we use to test our model. The flow simulations performed with a commercial flow simulator [13], used to compare the performance of our model, are described in Sec. IV. In Sec. V we compare the performance of our model with that of a commercial flow simulator. Conclusions and proposals for further developments are given in Sec. VI.

II. THEORY
The behavior of two or more fluids flowing through a porous medium is typically described with the multifluid extension of Darcy's law in combination with the relevant mass-balance equations for the different fluid fractions. The usual approach to modeling flow is to discretize these equations by a finite-volume approach so that the flow is then modeled on a grid where the rock and fluid properties are constant within each cell in the mesh. As noted in Sec. I, this can result in very-computationally-demanding simulations when one is attempting to model flow in geologically heterogeneous systems as millions of cells may be needed to describe the spatial permeability and porosity variations. In this paper we approximate the multifluid flow by first approximating the single-phase (tracer) flow through the geological system of interest (the sweep efficiency) using a growth model derived from that of Eden and then applying a first-order correction to approximate the effect of water displacing oil or non-aqueous phase liquid (NAPL) based on the analytical Buckley-Leverett solution [14] (the microscopic displacement efficiency). This approximation is reasonable provided that the viscosity ratio between oil and water is not too large. In this section we first describe our proposed growth model and then summarize how we estimate the first-order correction using Buckley-Leverett theory.

A. Growth model in geological systems
In its simplest form Eden's model grows clusters on a lattice beginning from a germ or seed at a single point or line of points on the lattice. The cluster grows by first occupying sites adjacent to the germ site(s) and then occupying sites adjacent to currently occupied sites. It is not possible to grow new clusters that are not connected to the original seed site(s), so this model can be seen to be at least partially analogous to Darcy flow (i.e., flow described on scales larger than those of individual pores) in porous media when the injected fluid is always connected.
The aim of our growth model is to estimate the area (two dimensions) or volume (three dimensions) of rock contacted by water when it flows through a rock initially filled with oil or NAPL with spatially varying permeability. As noted above, our model is based on the assumption that the two fluids have the same viscosity and that relative permeability effects are negligible (i.e., it describes tracer flow).
Following Eden's growth model, we describe the porous medium using a lattice or grid. The geologically heterogeneity is described by our assigning different permeabilities to different elements or cells in the mesh. This means that the permeability within each cell is constant but the pattern of permeability variation across the mesh approximates that seen in the geological deposit of interest [examples of such systems are given in Fig. 6]. This is consistent with the way that geological heterogeneity is modeled in finite-volume simulations of flow in oil reservoirs and contaminated aquifers. We use a regular Cartesian mesh, but the method can easily be extended to other types of structured or unstructured meshes. For single-fluid tracer flow, Darcy's law is written as where q is the flux (or Darcy velocity), ∇P is the fluidpressure gradient, μ is the fluid viscosity, and K is the permeability of the porous medium. We now apply Eq. (1) to the two situations represented in Fig. 1, which shows two systems with different injection-well geometries: a line source and a point source. In each case the injection well is drawn as a gray-shaded block connected to four cells, each with a different permeability K i , i = 1, 4. If we assume a constant pressure drop between the injection area and each block, we can see that if the injected flow rate is Q, the flow rate Q i into each cell will be proportional to its permeability: This is equivalent to calculating the overall resistance of a set of resistors in parallel in an electrical circuit (remembering that permeability can be treated as equivalent to electrical conductivity). We define our cluster-growth rules on the basis of Eq.
(2) such that the probability for water to move from an occupied cell to an unoccupied cell will depend on the FIG. 1. Two different well geometries represented on a regular Cartesian grid. Both geometries represent parallel-flow situations. It is possible to see from Eq. (1) that, assuming a constant pressure drop between the injection area and each block, the flow rate in each block will be proportional to its permeability as is shown by Eq. (2).
permeability of an unoccupied cell compared with that of its neighbors. Figure 2 illustrates how these growth laws are defined. The central cell acts as a source with a constant flow rate Q. According to Eq. (2) this flow will fill neighboring cells at a rate proportional to their permeabilities. We set a threshold to define when a cell is full. We can see that in this way the higher-permeability cell will be filled faster than the lower-permeability one. In our model we consider that when a cell is filled, it becomes a source itself (i.e., it becomes part of the cluster). Figure 3 summarizes the growth process. As noted above, the porous medium is represented by its permeability distribution on a lattice [ Fig. 3(a)]. The initial seed is located along the line x = 0, the position of the water injector. A counter C i (initially set to zero) is placed at each new growth site [ Fig. 3(b), step 1]. At the beginning of each time step, a growth site i is randomly chosen [ Fig. 3(b), step 2] and its counter C i is increased by an amount α proportional to its own permeability K i : where we can interpret the amount αK i as the amount of fluid flowing into the block each time it is selected and the counter C i is related to the fluid saturation of the block.
If this counter C i reaches a predetermined threshold C th , the site becomes fully occupied (filled) and its neighbors become active [ Fig. 3(b), step 5], otherwise a new active block is (randomly) chosen and the process continues. The The basics of our model. A full cell floods on its neighbors. The flow rate for each neighbor will be proportional to their permeability. In this way higher-permeability cells will be filled faster than lower-permeability ones. Steps followed for the model: A seed is set in the position of the flow source and all activated cells (adjacent to the seed) are set with a counter set to zero (step 1). One of the active cells is selected at random (step 1, red circle) and an amount equal to its permeability is added to its counter (step 2). If the counter value becomes bigger than C th , the cell becomes occupied (step 5) and cells adjacent to this cell are activated (step 6). Otherwise the process starts again (steps 3 and 4) by randomly selecting a new active-unit block. The breakthrough time is defined as the moment at which the cluster reaches the y = L side.
threshold value C th and α = 1 are constant throughout the system (i.e., they are global system parameters in the current implementation). In future work, they could be chosen so that the effects of gravity or viscosity ratio are also represented.
The process continues until some block belonging to the region defined as the sink (line y = L in this case) is occupied. We call this the "breakthrough time" and we measure the fraction P of cells in the system belonging to the connected cluster at this time. This value, P, is our measure of the system sweep efficiency.
It is easy to see that when it is applied to the homogeneous case, the model behaves as a growth model with noise reduction controlled by the parameter C th . This noise-reduction method, first proposed to deal with three-dimensional cluster growth, ensures that blocks that were occupied first are more likely to grow than those that were just occupied. The model can also be interpreted as an aging stochastic invasion percolation model: locally, the blocks with the highest permeabilities have the highest probability of becoming occupied, but the lowerpermeability blocks will also become occupied eventually.
In this work we set C th equal to 5 times the maximum permeability in the system to select a global system parameter. The value of this threshold does not affect significantly the value of the sweep efficiency obtained for the system but it does affect the speed of the calculation of the growth. A bigger threshold value will slow down the performance, while smoothing the growing front. On the other hand, a lower C th will result in a shorter execution time but a much noisier front. In any case, the impact of this selection can be easily evaluated, as the algorithm is much faster than a typical flow simulation. The only calculation at every time step is an addition and a comparison of the new counter value for the selected cell with the threshold C th . No probabilities are calculated and no steps are rejected. This makes the algorithm highly efficient. The algorithm scales with N , the number of grid blocks in the system.

B. Correction using Buckley-Leverett theory
When two immiscible fluids F 1 and F 2 (e.g., oil and water) flow simultaneously through the pore space, the two fluids interact with each other and the pore walls reduce the overall flow rate. This is represented by assigning a relative-permeability function to each fluid that varies with the volume fraction of that fluid in the porous medium [ Fig. 4(a)], and we write the flow rate for each  [15], we know that the intersection between the straight line starting at the the connate water saturation S wc and the tangent to the fractional-flow curve and the axis y = 1 gives the average water saturation in the system at the breakthrough time [14]. Subtracting S wc from this value, we obtain the pore volume fraction of water injected in the system at the breakthrough time. S or , residual oil saturation. fluid as In general, the pressure in each fluid is also different, giving rise to a capillary-pressure function, but this is usually important only on length scales on the order of meters. In this work we consider length scales on the order of hundreds of meters, so we assume that capillary pressure is negligible.
To analyze the displacement process in a onedimensional porous material we use Buckley-Leverett theory [14]. Consider Fig. 5(a): a quasi-one-dimensional homogeneous porous system of length L and porosity φ initially occupied by a fluid F 1 . The cross-section area a of the system is such that √ a L. A second fluid F 2 is injected into one end of the system. According to Buckley-Leverett theory [14], the typical saturation curve S F 2 of the injected fluid F 2 as it displaces fluid F 1 is shown in Fig. 5(b). This consists of a shock front (i.e., a saturation discontinuity between the phases), followed by a gradual transition between the shock-front saturation and the maximum possible saturation of fluid F 2 , S MF 2 = 1 − S F 1 R , where S F 1 R is the residual saturation of the fluid initially in place. The specific shape of this saturation curve is determined by the relative-permeability curves [ Fig. 4(a)] and viscosities of the two fluids.  5. (a) A porous material of length L and cross-section area a. A fluid F 2 , injected from side x = 0, displaces the fluid F 1 , initially in place. (b) Characteristic saturation curve for an injected fluid F 2 as it displaces fluid F 1 through the porous material. We can see the initial saturation of the fluid F 2 , the remaining saturation S F 1 R of the fluid F 1 initially in place, the shock front, representing the saturation discontinuity between the phases, and the gradual transition between the shock-front saturation and the maximum possible saturation for the injected fluid.
Using the Buckley-Leverett equation and applying the graphical technique of Welge [15], we can estimate the theoretical fractional pore volume W iT b of the injected fluid at the breakthrough time [ Fig. 4(b)]. The theoretical breakthrough time T bBL for the system can then be estimated as the necessary time to inject that amount of pore volume of fluid F 2 into the system: We can see that if we use a constant volumetric flow rate Q, the breakthrough time will be independent of the absolute effective permeability of the system. We can normalize the sweep efficiency estimated from the growth model by the breakthrough time calculated from Buckley-Leverett theory to estimate the breakthrough time that would be seen when a waterflood is performed in the geological model under consideration.

III. TEST CASES
Two different types of geological heterogeneity are studied to test the model. The first, simpler case is formed of two layers of contrasting permeability parallel to the flow [ Fig. 6(a)] to explore the impact of different permeability contrasts on the sweep efficiency. The permeability of the higher-permeability layer is kept constant at 100 mD and the permeability of the lower-permeability layer is varied between 100 mD (i.e. a constant-permeability map) and 10 −2 mD. The second type of heterogeneity is based on the socalled model 2 of the tenth SPE comparative solution project [16]. This is a synthetic but geologically realistic representation of heterogeneities found in a North Sea Brent sequence and has become a standard test vehicle for evaluating new methods for modeling flow in macroscopically, heterogeneous porous media. Although the SPE 10 model 2 is three-dimensional, we subdivide it into a set of 85 two-dimensional layers. This enables us to test the growth model more thoroughly. The first 35 layers (which represent the Tarbert formation) have a statistically homogeneous permeability distribution [ Fig. 6(b)], while the lower 50 layers (which represent the more heterogeneous Upper Ness sequence) show distinct, high-permeability channels set within a background of much lower permeability [ Fig. 6(c)]. Each layer is 670.7 m long, 365.8 m wide, and 0.61 m thick and has a constant porosity of 0.4. The SPE 10 model 2 enables us to test the growth model with a wide range of different permeability distributions, each one with a markedly different characteristic spatialcorrelation behavior. All geological models are discretized onto a (220, 60, 1) Cartesian grid. This grid is chosen to be consistent with the discretization of the SPE 10 model 2.

IV. FLOW SIMULATION
We perform full-physics flow simulations using a commercial black oil simulator [13] on all test cases to provide a benchmark against which to compare the predictions of our growth model. We simulate water injection into a reservoir containing oil. We use the relative-permeability curves given in Fig. 4(a). The oil and water viscosities are both 5 × 10 −4 Pa s, giving a viscosity ratio of 1. The water and oil densities are both 1033 kg/m 3 . It is assumed that the capillary pressure between phases is zero. An injection well is completed in all the cells along the left-hand face of each model and a production well is completed in all the cells along the right-hand face of the model. Constantrate injection and production with a voidage replacement ratio of 1 (as the fluids are incompressible) are applied. We record the breakthrough times for water obtained from the different simulations and normalize these by dividing them by the breakthrough time calculated using Buckley-Leverett theory [Eq. (5)] to properly compare them with the results from the growth model. Of course, we would actually multiply the growth-model sweep efficiency by the Buckley-Leverett breakthrough time if we were using it to replace conventional flow simulation.
In Fig. 4(b) we can see that the fractional pore volume injected at the breakthrough time is 0.606. For a volumetric flow rate of Q = 3 m 3 /day, this corresponds [Eq. (5)] to a 034003-5 breakthrough time of 11 886.5 days in our model reservoir if each one of the analysed systems were homogeneous.

V. RESULTS
In this section we compare the fraction P of the cluster obtained for the systems in Sec. III against the breakthrough time obtained with the multiphase flow simulation. Assuming a constant injection rate Q, the breakthrough time for a homogeneous system is given by Eq. (5). This time will correspond to a cluster fraction P = 1. All the breakthrough times obtained from the flow simulations are normalized with this theoretical value. It is important to mention here that although Eq. (5) is valid for homogeneous porosity, it is easily extended to the case of nonhomogeneous porosity. Figure 7 shows the behavior of the growing cluster and the behavior of the saturation of the injected fluid obtained from the commercial flow simulator. We can see that there is good agreement between the shapes of the patterns for the two cases and also that the growth model performs well for the different types of permeability distributions studied: layered [ Fig. 7(a)], statistically homogeneous type [ Fig. 7(b)], and highly spatially correlated [ Fig. 7(c)]. Figure 8 shows the correlation between the size P of the connected cluster obtained with our growth model and the normalized breakthrough time obtained with the commercial flow simulator [normalized with use of Eq. (5) and the system properties given in Sec. IV]. The straight line shows a proposed linear fit to the data. Figure 8(a) shows the result corresponding to the twolayer case for nine different K high /K low ratios. We can see that the linear fit gives a correlation very close to 1 (0.967 ± 0.007). The CPU time required for the growth model is on average 0.11 s per realization, while each fullphysics simulation takes on average 673 s per realization. Figure 8(b) shows the results from the 85 layers of the SPE 10 model 2. Blue triangles correspond to the statistically homogeneous cases [ Fig. 6(b)] and red circles correspond to the highly channelized ones [ Fig. 6(c)]. We can see very good agreement between the results obtained with our model and those obtained with the commercial simulator. The later breakthrough times and the bigger spread corresponding to the statistically homogeneous cases [ Fig. 8(b), blue triangles] are expected because flow has a wider selection of possible paths to take here. On the other hand, the channelized cases show a smaller spread [ Fig. 8(b), red circles] as the flow is much more focused here. For these cases the growth model requires on average 20 s per case, while the flow simulation requires 2060 s. All simulations are performed with a single processor.

VI. CONCLUSIONS
We present a growth model for the rapid estimation of sweep efficiency in geologically heterogeneous porous media. This model can be used in combination with the Monte Carlo method to assess the impact of permeability uncertainty on breakthrough time and sweep in hydrogeology and hydrocarbon recovery. It can also be used to estimate the effect of different upscaling methods. The model takes into account the local spatial correlation of permeability and statistically estimates patterns that the injected fluid could be more likely to follow.
The proposed model allows us to estimate the proportion of the total system that is swept by the injected fluid and the breakthrough time. We achieve precision similar to that of a commercial simulator in a fraction (namely, 1/1000-1/100 depending on the permeability pattern) of the time. This speedup would allow us to screen different scenarios quickly to characterize the effect of geological uncertainties on oil recovery or contaminant remediation in aquifers. The model does not require the pressure equation to be solved (unlike the case with models such as streamline models [17,18]) and does not need any extra scaling parameters. We show that it can handle different permeability distributions and spatial correlations provided these can be represented on a discrete mesh or grid. This is typically the case for both hydrogeological and petroleum applications. This shows that the algorithm presented is highly robust and versatile.
The algorithm is particularly efficient as the only calculation at each time step is the generation of a random number, and addition and a comparison of permeability values. We use a square lattice here but different grid geometries could be used. Although the sweep efficiency can be estimated for any boundary conditions (constant, rate, constant pressure, etc.) the breakthrough time can be estimated only by assuming a constant-flow-rate injection. It would be useful to address this issue in further work; namely, how to incorporate a timescale into the model when a constant pressure or another, different, boundary condition is used instead. It will potentially provide us with a powerful tool to estimate the well production decay after the breakthrough time. The model is demonstrated with two-dimensional geological models with isotopic permeability but can easily be applied to three-dimensional systems, provided gravitational effects are negligible. It should be possible to include the effects of gravity and permeability anisotropy by making α directionally dependent. Further work should also deal with how to incorporate different water-oil fluid properties, such as viscosity.