## versión On-line ISSN 1807-0302

### Comput. Appl. Math. vol.28 no.3 São Carlos  2009

#### http://dx.doi.org/10.1590/S1807-03022009000300003

The uncertainty effects of deformation bands on fluid flow

Michael PreshoI; Victor GintingI; Shaochang WoII

IDepartment of Mathematics, University of Wyoming, 1000 East University Avenue, Laramie, WY 82071
IIEnhanced Oil Recovery Institute, University of Wyoming, 1000 East University Avenue, Laramie, WY 82071. E-mails: mpresho@uwyo.edu / vginting@uwyo.edu / swo@uwyo.edu

ABSTRACT

Reservoir fractures and deformation bands are capable of affecting fluid flow and storage in a variety of ways. In terms of flow effects, we typically encounter an unchanged or increased permeability when considering flow parallel to a fracture, whereas we expect a noticeably reduced permeability when considering flow across (perpendicular to) a deformation band. For this paper, we refine our efforts and focus on the effects of deformation bands on multi-component porous medium flow. The main assumption is that the width of the band is a random (uncertain) variable that follows a certain statistical distribution from which a number of realizations can be generated. Monte Carlo simulations can then be performed to obtain statistical representations of the transport quantity in relation to the nature of uncertainty. As analytical expressions are available for this quantity of interest, we are able to compare them with the Monte Carlo results. Furthermore, we derive a stochastic perturbation model as an alternative to Monte Carlo simulations. Finally, a set of numerical examples is presented to illustrate the performance of these approaches.
Mathematical subject classification: 65C20.

Keywords: deformation band, pressure equation, saturation equation, Monte Carlo, uncertainty, stochastic perturbation expansion.

1 Introduction

Deformation bands typically represent sections of porous media with significantly reduced permeability and porosity resulting from cementation [8]. The decreased permeability inhibits fluid flow, contributes to a higher pressure drop, and as a result, deformation bands can play a role in reservoir fluid flow [8]. The study of deformation bands extends to various areas of academia and industry, and in this paper we address the topic within a mathematical framework [1, 8, 17, 18]. For our particular petroleum industry application, we are interested in quantifying the effects of deformation bands on oil production. This can be done by assigning uncertainty to various parameters that describe a band. The effects of deformation band permeability, orientation, and width (aperture) have been studied in a non-mathematical framework [8]. As a result, we assume that the above parameters can represent random variable candidates in the initial problem formulation.

In this study, we treat the width of a deformation band as the random variable of interest. This is a natural choice since variation in the width of a band is often guaranteed in a subsurface reservoir [8]. In other words, knowledge that deformation bands exist in a reservoir does not offer much information about the width(s) of the bands. We also point out that wider bands progressively inhibit flow (due to the reduced permeability), and thus the width is a parameter that can most directly affect oil production in a porous medium [8]. With width as the random variable, we make deterministic assumptions on the permeability and location of a vertically oriented deformation band in an idealized porous medium. These initial assumptions offer a foundation in describing the physical model. The multi-component flow under consideration is modeled by a hyperbolic partial differential equation coupled with Darcy's Law [2, 3, 4, 16]. Within this setting, the velocity and saturation become functions of a random variable due to their dependence on the deformation band's width. In turn, the oil production is also random.

By appending a source of uncertainty to the width of the band, we can analyze the oil production in a variety of ways. To address the uncertainty of the width we initially assume the parameter is governed by a given probability density function. By considering idealized models that offer the luxury of analytical solutions we can use the probability density function to compute exact statistical mean values. These exact values give a benchmark of comparison with numerical methods that assess the effects of uncertainty. One such method is the Monte Carlo method. This method involves generating a number of realizations of a random variable and then averaging over the realizations to obtain effective solutions [9]. In the scope of our problem, we choose a probability density function, generate a number of realizations for the width of a deformation band, and compute the statistical mean production curves. Monte Carlo is a widely used tool and much of today's computing power is devoted to calculating Monte Carlo solutions [6]. However, one main disadvantage is that Monte Carlo is a very expensive method to implement, particularly if a large number of realizations is required. It is this disadvantage that leads us to explore more efficient methods of computing expectations.

In order to address the costly aspect of Monte Carlo, we introduce a stochastic perturbation model. The derivation of the model hinges on the assumption that the saturation and velocity can be expressed as first order stochastic perturbations (see, for example [13, 14, 19]). Using this assumption we can then reintroduce the variables into the original equations to obtain a governing equation for the saturation statistical mean. In deriving the model we make two main assumptions. First, we assume that a first spatial derivative does not change significantly along a characteristic. Second, we assume that any third order stochastic terms can be neglected. The resulting model is, in general, more efficient than Monte Carlo and offers solutions that very closely match Monte Carlo for low deviations. We note that this approach has been previously introduced by Glimm and Sharp, and Zhang [11, 21] (see also [20]). Within the context of upscal-ing in heterogeneous flow a similar approach has been used in [5, 7, 10, 15]. However, for higher deviations we encounter solutions that stray slightly from Monte Carlo. This is due to the fact that we neglect a third order stochastic term in the model derivation. We classify the latter results as a limitation of the model, however, we remark that the model still performs reasonably well regardless of some discrepancies found for higher deviations.

In §2 we introduce the physical transport model and coupled pressure equation for multi-component flow. Under the assumption that the width of the deformation band is random we solve the respective set of equations analytically. By doing so, we eliminate any errors that may arise from numerical approximation and focus solely on the effects of uncertainty. In §3 we use the analytical solutions to calculate the Monte Carlo production curves using both Uniform and Gaussian distributions [12]. As the exact statistical mean is available in these cases, we can make a comparison with the Monte Carlo results. In §4 we derive a stochastic perturbation model which hinges on expressing the saturation and velocity as first order stochastic perturbations. In §5 we solve the model numerically and compare the results with the Monte Carlo solutions from §3. In addition, we perform a probability density function comparison. Finally, we offer some concluding remarks in §6.

2 Model equations

In this section we begin by introducing the hyperbolic saturation equation that models a general two-component system. Letting S denote the saturation of a displacing fluid, we consider

where x [0, L], t [0, ),and υ(t) = is the flux obtained through the pressure equation

which is coupled to (1) through the viscosity

The parameter is the viscosity ratio, μo denotes the viscosity of the oil component and μd denotes the viscosity of a displacing fluid. In the general two-component flow, we assume that M < 1. This models a polymer flood situation where the polymer (displacing fluid) component is more viscous than the oil component. A single phase flow is a special case in which, M = 1, i.e., there is only one fluid present.

Within this framework, we use βto denote the random width of the deformation band and assume a piecewise constant permeability structure

See Figure 1. We are particularly interested in the situation wherek1 = k3 =kfandk2 =kb < kf. Here kfrepresents the field permeability, and kb represents the permeability within a deformation band.

We point out that β, is now explicitly present in the problem, and combining all these equations ensures that the resulting saturation solution is random. We can describe the general idea in a simple flow chart: β k(x; β)v(t;3)S(x, t; β).

In solving (1) analytically we ultimately want to derive an expression for the location of the saturation front. We denote this location by ξ(t; β).Upon successfully finding this quantity, the saturation solution is

Furthermore, we eventually want to analyze the uncertainty effects of the production curve F (t; β)defined as

The production curve essentially shows us the saturation ofoil at the right boundary, L, ofour interval. See Figure 2 for an example ofa production curve plotted against time.

Generally speaking, to solve for ξ(t; β),we must first derive the flux. Then we set . From here we arrive at an equation for ξwhich can easily be solved. Initially, it is important to note that we must consider a number of cases for the flux and front location. See Figure 3. In particular, we have three separate flux values

where x0 =0, x1 =, x2 =, and x3 = L depending on where ξ lies in the spatial interval. In addition, when the porous medium is fully saturated,

i.e., the front has passed the right boundary. Treating each case separately, using (2) and (4), and ensuring continuity at ξwe obtain

where Δp= p0 — pL, and mi,ni, and qi are deterministic values depending on the parameters kf, kb, L, and M. As a representative example,

Furthermore,

is the fully saturated flux value, where ns and qs are similarly deterministic. Using (9) we can now set . Essentially, we have a first order, separable differential equation for each case. Treating each case separately and ensuring continuity we arrive at a quadratic equation for ξi for i = 1, 2, 3. In addition, we can solve for the βdependent times when the front values shift:

See Figure 4 for a plot of the front locations and time values. Of particular interest, ξ3 takes the form

Then (13) is appropriately used to evaluate S(L, t) from (1). We remark that we are considering a model that yields a random fluid velocity. In turn, (1) gives a random, analytical saturation solution. With this solution we can now implement the Monte Carlo method and compute statistical mean values of production.

3 Monte Carlo results

The Monte Carlo simulation is implemented for two-component and single phase flows. We generate N = 10, 000 positive realizations of β from Gaussian and Uniform distribution, with fixed bñ= 10 and various standard deviations, σb. For each realization βi we compute ξ3i (t; β), i = 1, ..., N from (13), and subsequently a production curve Fi (t) from (1). We then average the values at each time level and plot the resulting statistical mean production curve given by

The deterministic data are kf = 500, kb = 100, L = 100, p0 = 1000, pL = 0, μo= 2.7, and μd = 3 (for two-component flow) and μo = μd = 1 (for single phase flow).

For the problem we consider, convergence of the Monte Carlo simulation can be accessed by way of comparison to the exact production curve statistical mean expressed as

where fβ(β) is the respective probability density function of the Gaussian or Uniform distribution [12]. This comparison is presented in Figure 5, showing for single phase flow and two-component flow. We chose σb= 10 as a middle ground for comparison. The two results are indistinguishable from each other which indicates the convergence of Monte Carlo for theN= 10, 000 realizations.

Next, we present in Figure 6 the behavior of the production curve for various standard deviation σb. All plots indicate that for increasing σb's we see a more predominant deviation from the deterministic solution (σb=0) near the breakthrough times. For realizations from Gaussian distribution, we see solutions that predictably deviate smoothly in time. However, we see sharper drops to the left of breakthrough time for higher deviations. We generally expect earlier breakthrough when assuming more uncertainty. In the case of the realizations from the Uniform distribution, we see solutions that sharply deviate from the deterministic solution. Higher deviations lead to predictably increased variation after breakthrough and earlier breakthrough times. In general, we see similar pattern for higher deviations only with sharper curves.

We note that for higher deviations, negative realizations are more likely occur [12]. This is physically unrealistic for the application under consideration. Thus, in the process we exclude all these negative values in the Monte Carlo simulation. Consequently, this leads to a situation where the Monte Carlo curves are non-symmetric with respect to the σb =0 curves.

4 Stochastic perturbation expansion

Thus far we have computed statistical mean of the production curves using Monte Carlo simulation and gave comparison with their exact counterparts. We reiterate the fact that Monte Carlo is, in general, an expensive computational approach. This is particularly the case when a large number of realizations is desired. To address this issue we now consider a separate, less costly method of computing statistical values in relation to our problems. In particular, we will derive an equation governing the statistical mean of the saturation S using the notion of Stochastic perturbation expansion. Recall the fact that v(t) and S(x, t) are random functions due to their dependence on β. To simplify notations, we will no use

Stochastic perturbation assumption allows us to express these variables as

where S'(x, t) and v'(t) are stochastic fluctuation terms. Substitution of (15) into (1) and appropriately rearranging the resulting equation yield

By taking the expected value of (16) and using 0, we get

Here we have obtained a governing equation for statistical mean of the saturation. Our next task is to model the second order effect To accomplish this we now work on the characteristics Using the notion of a total derivative, (16) becomes

Integration along the characteristic, with the assumption that does not significantly change along the characteristics yield

Multiplying by υ'(t), taking expectation, and neglecting the higher order stochastic term give

and therefore,

Substitution of (19) into (17) yields

where

For single phase flow, the flux is independent of time and thus The above differential equation is completed by imposing a natural boundary condition in addition to the existing initial and boundary conditions from the original problem.

We note that more general equations were derived in [11] and [21] for single phase flow. In particular, Glimm and Sharp (see [11]) and Zhang (see [21]) obtained related results in the context of a multi-length-scale, random permeability field. For two phase flow, we refer the reader to [20].

5 Comparison of approaches

This section is devoted for comparison between the Monte Carlo simulation and the Stochastic Perturbation Expansion described in the previous section. To solve (20) we use a Backward Euler, centered finite difference scheme where

The equation is solved on a very fine grid and very small time steps in order to minimize any sources of error from numerical discretization. In particular, we want to focus solely on the effects of uncertainty. We point out that the first step in solving the model is the computation of and a(t).

Comparison for single phase flow is presented in Figure 7. We use a relatively low deviation, σb =2, and a relatively high deviation, σb =10. In addition, the Gaussian and Uniform distributions are used for modeling the random behavior of β.

We point out that for σb =2 in the Gaussian case, the numerical solution and the Monte Carlo solution coincide nearly exactly. Thus, the model accurately portrays the Monte Carlo solution for a low deviation. This also indicates that the assumption that does not significantly change along characteristics in §4 is reasonable. For σb =10 we see an accurate portrayal of the breakthrough time, yet there is a noticable smoothing at the drop. This is expected as the model is parabolic. We recall that in §4 we also made the assumption that a third order stochastic term could be neglected. For σb =10 we simply see the natural breakdown of this assumption.

We now shift our attention to the Uniform case. In this case, we encounter similar behavior of the model solutions (with respect to the Gaussian model solutions). The main difference is that we see pronounced smoothing as compared to the Uniform Monte Carlo curves. Again, this is expected since the model is parabolic, whereas the Monte Carlo curves are inherently sharper in the Uniform case. Similarly, we see an accurate portrayal of breakthrough time even for a higher deviation, yet a discrepancy after that.

The results for the model applied to the two-component case are shown in Figure 8. We use a relatively low deviation, σb =2, and a relatively high deviation, σb =10 for the Gaussian disbribution, while for the Uniform distribution we use σb =3 and σb =10.

The stochastic perturbation model agrees with the Monte Carlo solutions for a low deviation, while similar behavior as in single phase flow occurs for a higher deviation. We encounter similar behavior of the model solutions for the Uniform distribution (see Fig. 8). The main difference is that we see pronounced smoothing as compared to the Monte Carlo curves. Again, this is not a surprise since the model is parabolic, whereas the Monte Carlo curves are inherently sharper in the Uniform case.

Next we offer a comparison between the Gaussian and Uniform distribution results. At this point, we have treated each case separately in computing the Monte Carlo and numerical solutions. However, we think it is beneficial to offer a side by side comparison. For these comparisons we only consider the two-component case, although the same conclusions hold in the single phase case.

The left side of Figure 9 compares the respective Monte Carlo solutions for the arbitrary case when σb =10. As mentioned earlier, the Gaussian distribution gives inherently smoother Monte Carlo solutions. Nevertheless, we see that the breathrough time is identical regardless of what distribution we use.

We also compare the solutions of (20) as solved with both the Gaussian and Uniform distributions (right side of Figure 9). Interestingly, the numerical solutions are quite similar to one another. This indicates that the model is not very sensitive to the probability density function we use to calculate the required coefficients. The curves are identical near the breakthrough time, and differ ever so slightly after breakthrough.

6 Conclusion

In this paper we introduce a stochastic perturbation model to quantify the effect of uncertainty of the deformation band width to oil production. By expressing our saturation and velocity quantities in terms of first order stochastic perturbations we are able to obtain from the original equations a parabolic equation modeling the statistical mean values of saturation. For low deviations we encounter solutions that almost exactly coincide with the Monte Carlo results. For higher deviations we see discrepancies with the Monte Carlo results that can be attributed to the assumptions in the model derivation. We ultimately obtain a more efficient method for computing expected saturation values and we verify the accuracy by first computing Monte Carlo production curves (with a variety of deviations) for the single phase and two-component saturation equations. In addition, we conclude that the model is not sensitive with respect to the probability density functions of interest (Gaussian or Uniform).

REFERENCES

[1] A. Aydin, Small Faults Formed as Deformation Bands in Sandstone. Pageoph, 116 (1978), 913-930.         [ Links ]

[2] J. Bear, Dynamics of Fluids in Porous Media. Courier Dover Publications (1988).         [ Links ]

[3] L.P. Dake, Fundamentals of Reservoir Engineering. Elsevier (1978).         [ Links ]

[4] H. Darcy, Les Fontaines Publiques de la Ville de Dijon. Dalmont, Paris (1856).         [ Links ]

[5] Y.R. Efendiev, L.J. Durlofsky and S.H. Lee, Modeling of subgrid effects in coarse scale simulations of transport in heterogeneous porous media. Water Resour. Res., 36 (2000), 2031-2041.         [ Links ]

[6] D. EstepandD.Neckels, Fast methods for determining the evolution of uncertain parameters in reaction-diffusion equations. Comput. Method Appl. Mech. Eng., 196 (2007), 3967-3979.         [ Links ]

[7] R. Ewing, Y. Efendiev, V. Ginting and H. Wang, Upscaling of Transport Equations for Multiphase and Multicomponent Flows. Lect. Notes Comput. Sc. Eng., 60 (2008), 193-200.         [ Links ]

[8] H. Fossen and A. Bale, Deformation bands and their influence on fluid flow. AAPG Bull., 91(12)(2007), 1685-1700.         [ Links ]

[9] J.E. Gentle, Random Number Generation and Monte Carlo Methods. Springer Verlag, (2003).         [ Links ]

[10] V. Ginting, R. Ewing, Y. Efendiev and R. Lazarov, Upscaled modeling in multiphase flow applications. Comput. Appl. Math., 23(2-3) (2004), 213-233.         [ Links ]

[11] J. Glimm and D. Sharp, A Random Field Model for Anomalous Diffusion in Heterogeneous Porous Media. J. Stat. Phys., 62(1-2) (1991), 415-424.         [ Links ]

[13] M. Kaminski, Application of the generalized perturbation-based stochastic boundary element method to the elastostatics. Eng. Anal. Bound. Elem., 31(6) (2007), 514-527.         [ Links ]

[14] M. Kammski, Generalized Stochastic Perturbation Technique in Engineering Computations. J. Eng. Appl. Sci., 3(3) (2008), 246-260.         [ Links ]

[15] P. Langlo and M.S. Espedal, Macrodispersion for two-phase, immisible flow in porous media. Adv. Water Resour., 17 (1994), 297-316.         [ Links ]

[16] P.D. Lax, Hyperbolic Systes of Conservation Laws and the Matheatical Theory of Shock Waves. SIAM, Philadelphia (1973).         [ Links ]

[17] T. Manzocchi, P.S. Ringrose and J.R. Underhill, Flow through fault systems in highporosity sandstones. Struct. Geol. Reservoir Characterization, 127 (1998), 65-82.         [ Links ]

[18] K.R. Sternlof, M. Karimi-Fard, D.D. Pollard and L.J. Durlofsky, Flow and transport effects of compaction bands in sandstone at scales relevant to aquifer and reservoir manageent. Water Resour. Res., 42 (2006), W07425.         [ Links ]

[19] D. Tartakovsky and A. Guadagnini, Prior mapping for nonlinear flows in random environments. Phys. Rev. E., 64 (2001), 03530(R).         [ Links ]

[20] D. Zhang, Stochastic Methods for Flow in Porous Media. Coping With Uncertainties. Academic Press, San Diego, CA (2002).         [ Links ]

[21]Q. Zhang, A Multi-Length-Scale Theory of the Anoalous Mixing-Length rowth for Tracer Flow in Heterogeneous Porous Media. J. Stat. Phys., 66(1-2) (1992), 485-501.         [ Links ]