A Mathematical Journal Vol. 6, No 4, (53-72). December 2004. On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs1 Frederico Furtado Department of Mathematics, University of Wyoming Laramie, Wyoming 82071-3036, U.S.A. furtado@uwyo.edu Felipe Pereira Instituto Politécnico, Universidade do Estado do Rio de Janeiro Rua Alberto Rangel, s/n, Nova Friburgo, RJ, 28601-970, Brazil pereira@iprj.uerj.br ABSTRACT The basic goal of scale up procedures is to simulate on coarse grids, with modified equations, multiphase reservoir flow and transport problems defined on fine grids. Scaling up is in general difficult. Difficulties arise both from the highly nonlinear fluid-fluid interactions in the flow of multiphase fluid mixtures and from the complex interactions between heterogeneities and nonlinearities. Our recent results show that several different flow regimes occur, depending on the relative strengths of flow nonlinearity and medium heterogeneity, as well as on the spatial structure of such heterogeneity. In the present study such results are used to investigate the applicability of a simplifying assumption made in several studies in the development of theories for the scale up problem for immiscible (water-oil) displacement. Our results indi- cate that such assumption is not appropriate to describe a typical mixing regime encountered in petroleum reservoirs, namely, the nonlinear unstable regime, in which nonlinearities in the governing equations dominate the fluid mixing process. Key words and phrases: Porous media, scale up, multiphase flow, heterogeneity, fractals, scaling laws 1This work was supported in part by the National Science Foundation (NSF) under Grant INT- 0104529 and ANP/CNPq under the “Plano Nacional de Ciência e Tecnologia do Setor Petróleo e Gás Natural (CTPETRO)”. 54 Frederico Furtado and Felipe Pereira 6, 4(2004) 1 Introduction Stochastic partial differential equations arise in theories which model complex (ran- dom) physical systems at small (microscopic) length scales, while their goal is to provide a description of the behavior of the system at large (macroscopic) length scales [28]. Thus, an essential step in these theories is the derivation of effective (or averaged) equations along with scaling laws governing their solutions, capturing essential features, at the macroscopic level, of the random microscopic structure. In mathematical studies of the transport of pollutants in groundwater and of oil recovery processes one is faced with the situation described above, where a system of stochastic partial differential equations models the two-phase flow in a porous medium. The system is composed of two equations, a transport equation for the saturation (the relative volume of one of the two fluids) coupled to an equation for the velocity field, which is given by Darcy’s Law and the incompressibility condition for the flow. The randomness enters the problem through the unknown properties of the rocks, especially the permeability tensor. A question of a major importance in applications concerns the growth of the region where mixture of two fluids occurs, which is produced by multi length scale rock heterogeneity as well as by fluid instabilities. We focus our attention on the length of the mixing region (the mixing length). We are specially interested in the study of anomalous (non Fickean) diffusive mixing. In the case of oil recovery, anomalous diffusion leads to early breakthrough of the water drive and reduced recovery, while in groundwater flow anomalous diffusion leads to the rapid growth of contaminant plumes. Statistical models have been proposed to describe rock heterogeneity properties which are incompletely known, especially at length scales below the interwell spac- ing. We consider fractal (or self similar) models for geological variability, which were introduced in earlier theoretical work [29, 4, 36, 27]. The use of fractal models is a basic advance over earlier modeling based on fixed length scale heterogeneities. The self-similar hypothesis is very powerful. It allows reconstruction of heterogeneity on all length scales from measurements on a small range of length scales, e.g., labora- tory and pilot-scale field studies, and provides improved agreement with the observed scale dependent dispersivities in the field data and with the observed geological het- erogeneities on all length scales. There is, however, no known law of geology which relates variability on distinct length scales. We view self similarity as a convenient tool to create simulated data which encompasses many significant length scales. As geological data appears to be nonuniversal and slowly varying as a func- tion of length scale, the multi fractal (rather than pure fractal) approach [46] has been introduced. The basic assumption of the multi fractal approach is that data varies smoothly in log-log variables, having well defined tangents with nonuniversal (scale dependent) slopes. Such an assumption provides a very general description of variability, in contrast to the self similarity assumption, and may provide better mod- els for specific reservoirs or aquifers. In order to investigate the effect of rock heterogeneity alone, the authors and collaborators have considered the passive tracer flow problem (linear transport). The 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 55 determination of the dynamics governing the growth of the mixing region for such model is an open problem as far as full rigorous proofs are concerned. However, a deep insight into this problem and evidence of the existence of asymptotic scaling laws governing the fluid mixing can be achieved through perturbative methods and high resolution computations, used to solve the stochastic system of differential equations. A characterization of the long time asymptotic growth of the mixing length was obtained in Refs. [27, 46] within distinct analytic approximations at the level of perturbation theory. These theories, developed in the Eulerian picture, give the time asymptotic growth rate of the mixing length as a function of the distance asymptotic exponent of the permeability field. The finite time transient effect was also studied in Ref. [46]. These theoretical predictions are restricted to the regime of small fluctuations of either porous media heterogeneity [27] or velocity field [46]. High resolution numerical simulations are not limited to the small fluctuations regime. In the absence of an exact solution for the problem of determining the time- dependent behavior of the mixing length, numerical simulations provide the only available tool in placing limits on the validity of the perturbative analysis and in the prediction of complex fluid flow behavior outside the regime governed by perturbation theory. The tracer flow problem was carefully investigated numerically in Refs. [38, 25, 26, 18, 19, 20] and the predictions of perturbation theory were fully validated. To add more references to these, in the context of linear flows, we refer the reader to Ref. [2] for a renormalization group approach, to Ref. [44] for a perturbation expansion analysis, and to Ref. [1] for a homogenization treatment. See also Refs. [13, 12, 10, 11, 8, 23, 35, 47, 48] and the references cited therein. Although still an area of active research, the body of work described above lead to a very good understanding about the scale up problem for linear transport problems in heterogeneous formations. For nonlinear problems the study of the scale up problem is considerably more difficult, due to the complex interactions between heterogeneities and nonlinearities. Consequently, nonlinear flows have not been so extensively investigated. In most existing approaches for scaling up multiphase (nonlinear) flow systems only the first order (hyperbolic) transport terms are modified. This is known as hyperbolic renormalization, or “pseudoization” of the laboratory scale relative per- meability curves. For a review of some standard procedures of this type we refer the reader to Refs. [5, 7, 16]. Effective procedures have been proposed recently in the literature [40, 41, 42]. Alternative scale-up methods, not based on a hyperbolic renormalization, can be found in Refs. [17, 32, 9, 14, 45, 30, 31]. Recent results of the authors (see Ref. [22]) show that several different flow regimes occur, depending on the relative strengths of flow nonlinearity and medium heterogeneity, as well as on the spatial structure of such heterogeneity. See also Refs. [3, 24, 33, 43] for related work. It is then plausible to expect that different scale-up methods (different types of coarse-scale models), tailored to the specifics of the distinct flow regimes, might provide a better coarse-grained description of multi- phase flow in heterogeneous petroleum reservoirs. Thus we believe that an improved understanding of the interplay between heterogeneity and nonlinearity, as we have 56 Frederico Furtado and Felipe Pereira 6, 4(2004) pursued, will furnish essential insight about existing techniques for scaling up and their limits. In particular, in this paper we investigate the validity of a simplifying assumption made in the theoretical developments reported in Refs. [30, 31, 45, 32, 14] through high resolution numerical simulations. In these studies it is assumed that the transport equation for the saturation is not coupled to the equation for the velocity field. Our results indicate that such assumption is not appropriate to describe the flow regime dominated by nonlinearities (or, equivalently, the regime characterized by weak heterogeneities). This paper is organized as follows. The stochastic model for two-phase immiscible displacement considered in this work is described in Section 2. In Section 3 we discuss our strategy for solving numerically the model introduced in Section 2. The method we have developed to perform a quantitative analysis of the macroscopic water-oil mixing process is presented in Section 4. The recent effort of the authors in identifying a new family of mixing regimes which were investigated through a high resolution numerical study of the problem at hand is reviewed in Section 5. We then turn to the main result of this work in Section 6, with a discussion of the applicability of recently proposed theories for the scale up problem for immiscible (water-oil) displacement. Our conclusions along with a discussion of important open problems related to the material discussed here will appear in Section 7. 2 Stochastic Modeling of Two-Phase Flow In this paper we consider the scale up problem for two-phase, immiscible, incompress- ible displacement in petroleum reservois. We neglect the effects of gravity, compress- ibility and capillarity and set the porosity equal to a constant. The two phases will be referred to as water and oil, and indicated by the subscripts w and o, respectively. The governing system of equations in the laboratory scale (see Ref. [37]) can be described as follows. Darcy’s law for each phase takes the form: ui = − kri(s) μi k(�x)∇p, where ui is the phase velocity, k is the absolute permeability, kri is the relative permeability and μi is the viscosity, each with respect to phase i, s is the water saturation, and p is the pressure. Darcy’s law combined with the conservation laws for the phases can be expressed in the familiar form of “pressure” and “saturation” equations: u = −λ(s)k∇p, ∇ · u = 0, (1) ∂s ∂t + ∇ · (f (s)u) = 0. (2) (The constant porosity has been scaled out by a change of the time variable.) Here, λ is the total mobility, f is the fractional flow of water, and u is the total velocity. 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 57 These parameters are given as λ(s) = krw(s) μw + kro(s) μo , f (s) = krw(s)/μw λ(s) , u = uw + uo. We consider Eqs. (1)–(2) in a two-dimensional rectangle Ω = (0, Lx) × (0, Ly), with the boundary conditions u · n = −q, on x = 0, p = 0, on x = Lx, u · n = 0, on y = 0, Ly, (3) where n is the outward-pointing normal vector to ∂Ω, and a uniform initial condition s(�x, 0) = s0. (4) Water is injected uniformly (at a constant rate q) through the left vertical boundary (x = 0) of Ω, no flow conditions are imposed along the horizontal boundaries (y = 0, Ly), and fluid is produced from a well kept at constant (zero) pressure at the right vertical boundary (x = Lx). The relative permeabilities are assumed to be: kro(s) = (1 − (1 − sro)−1s)2, krω(s) = (1 − srw)−2(s − srw)2, where sro and srw are is the residual oil and water saturations, respectively. In our studies we consider scalar, log-normal permeability fields, k(�x) = k0e ρξ(�x), (5) where ξ is a stationary Gaussian random field, characterized by its mean 〈ξ〉 = 0 (angle brackets denote ensemble averaging) and its covariance function C(�x, �y) =< ξ(�x)ξ(�y) > . The mean 〈k〉 and variance σ2k of the log-normal field k are set by the coefficients k0 and ρ. Changing ρ varies the coefficient of variation, CVk ≡ σk < k > , (6) of the permeability field. We use the coefficient of variation as a dimensionless measure of the heterogeneity of the permeability field. We take the field ξ to be isotropic and, to introduce variability over all length scales, fractal, or self-similar. Thus its covariance function is given by a power law: C(�x, �y) = |�x − �y|β, β < 0. (7) The scaling exponent β, known as the Hurst exponent, controls the degree of mul- tiscale heterogeneity: As it increases, the heterogeneities concentrated in the larger length scales are emphasized and the field becomes more regular (locally). 58 Frederico Furtado and Felipe Pereira 6, 4(2004) 3 Numerical Approximation of Two-Phase Flow 3.1 Computational solutions for flow and transport An operator splitting technique is employed. The saturation equation (2) and the pressure equation (1) are solved sequentially with distinct time steps. Tipically, for computational efficiency, larger time steps are used for the pressure calculation. Figure 1: Saturation surface plots displayed as a function of the viscosity ratio M for flow in a β = −0.5 permeability field with coefficient of variation CVk = 0.99. The values of M , from top to bottom, are M = 1, 2.657, 5. A detailed description of the numerical method that we employ for the solution of Eqs. (2)–(1) is given in [49] (see also [15] for an earlier version of this procedure). This 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 59 method has proved to be computationally efficient in producing accurate numerical solutions for two-phase flow problems. A second order accurate finite difference scheme (see Ref. [34]) is used for the dis- cretization of the saturation equation (2). This method can accurately resolve sharp fronts in the fluid saturations without introducing spurious oscillations or excessive numerical diffusion. For the global pressure solution, we use a (locally conservative) hybridized mixed finite element discretization equivalent to cell-centered finite differences [15], which effectively treats the rapidly changing permeabilities that arise from stochastic geology and produces accurate velocity fields. Important features of the simulation code described in [49] include mass conserva- tion of fluid phases, absence of grid orientation bias and reduced numerical diffusion. We illustrate the numerical solutions obtained with this simulation code in Figure 1. Saturation surface plots displayed as a function of the viscosity ratio M ≡ μo/μw for flow in a β = −0.5 permeability field with coefficient of variation CVk = 0.99. The values of M , from top to bottom, are M = 1, 2.657, 5. 3.2 Computational generation of fractal fields The (log) permeability field (a correlated Gaussian random process) is the image under a convolution mapping of a white noise (uncorrelated Gaussian) field. In more detail, let η = η(�x) be the white noise random field. The correlation function of η is < η(�x)η(�y) >= δ(�x − �y). (8) Let f ∗ g denote the convolution product of f and g, ( f ∗ g ) (�x) = ∫ f (�x − �y)g(�y)d�y, (9) and let f ∨(�x) = f (−�x). The permeability field fluctuation is given by ξ = f ∗ η, (10) where f is the convolution kernel to be specified below. By direct computation, the correlation function (covariance) of the field ξ is given by < ξ(�x)ξ(�y) >= ∫ f (�x − �z)f (�y − �z)d�z = ( f ∗ f ∨ ) (�x − �y). (11) To get the desired asymptotic behavior given by (7) the convolution kernel is taken to be f (�x) = c1 (c2 + |�x|)τ · c3 (c3 + |�x|)d/2 , (12) where τ = (d − β)/2, d is the number of space dimensions (d = 2 in our case), and c1, c2 and c3 are constants which set the overall strength of the correlation function, regularize its behavior near the origin, and regularize its behavior near infinity. For 60 Frederico Furtado and Felipe Pereira 6, 4(2004) τ < d (and x �= y), the convolution integral (11) is convergent at z = x and z = y, and we can take the limit c2 → 0 directly in (10) and (11). Similarly for d < 2τ , (11) is convergent at infinity, and we can take the limit c3 → ∞ in (10) and (11). Since we are interested in the range β < 0, we have d < 2τ and c3 = ∞ in all cases, while τ < d corresponds to the restriction β > −d. Either β > −d or c2 > 0 is required for the covariance to be locally integrable, and because the covariance is positive, this condition is also required for the covariance to be a distribution, as is necessary in the framework of generalized random fields [39]. Thus we require β > −d or c2 > 0. We now discuss the spatial discretization of the procedure just explained. The basic white noise Gaussian field η is now defined on a discrete level, i.e. is taken to be piecewise constant on the mesh blocks of a finite lattice. We distinguish two discretization errors, namely the local ones due to the finite mesh size of the η field, and the long range ones due to the finite extent of the η field computational grid. Our first results were to understand and then to minimize these errors. To do this, we studied the covariance of the discrete random field ξ = f ∗ η, which is just the discrete convolution product f ∗ f ∨, as compared to (7). Figure 2: Two examples of realizations of permeability fields used in the fluid flow simula- tions. The top picture refers to an uncorrelated Gaussian field (β = −∞) and the bottom one to a fractal field with β = −0.5. Cutoff values are 0.3 md for black and 4.0 md for white. We restrict our attention to β ∈ (−2, 0). In this case, the convolution (10) is not singular, and we can set c2 = 0, c3 = ∞ in (12). However we found that nearly singular integrals cause slow convergence, a fact which is further complicated by the 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 61 requirement that to analyze scaling laws over multiple length scales, very small effects are being observed. To accelerate convergence, we use nonuniform grids (local mesh refinement) to allow greater resolution of the local and long distance behavior of the convolution defining ξ and its covariance. To overcome the local and grid boundary errors, we introduce variable sized grids, allowing for finer local resolution and greater separation from the boundary. The spatial grids depend on β. As β → 0, long range blocks (large blocks) are used, while as β → −2, local refinement is important. We refer the reader to [38] for more details. The computed correlated field is illustrated in Figure 2 where two realizations of random permeability fields are displayed, corresponding to β = −∞ (uncorrelated Gaussian), top picture, and β = −0.5 (correlated fractal), bottom picture. 4 Quantitative Analysis of the Mixing Process A quantitative analysis of the mixing process is afforded by analysis of the growth rate of the mixing region as a function of time or, equivalently, travel distance. For this pur- pose, we introduce a time dependent length scale—the mixing length—characteristic of the extent of the mixing region and study its growth as a function of time. Figure 3: Saturation level curves for two-phase flow with viscosity ratio M = 5 in a β = −∞ permeability field (ragged contours). Superposed is the planar saturation profile for the corresponding homogeneous flow. 62 Frederico Furtado and Felipe Pereira 6, 4(2004) Since the bulk incompressible flow is driven by a constant injection rate, the fluc- tuations in u around its mean 〈u〉 = (q, 0) result in the spreading of the heterogeneous saturation fronts around their mean position, which is specified by the location of the planar homogeneous front. See Figure 3. An estimate of the spreading in any realization Kr of the permeability field is provided by �r = �r(t): �r(t) ≡ 1 (s− − s+) ∫ Lx 0 |s̄r(x, t) − SH (x, t)|dx. (13) Here, s̄r is the average in the direction transverse to the mean flow of the saturation solution corresponding to K = Kr; SH is the homogeneous saturation solution cor- responding to the constant permeability KH ≡ 〈K(�x)〉; s− (respectively s+) is the saturation value immediately behind (resp. ahead) the saturation front in SH . In the stochastic context, �r is a random variable, specified in terms of its statistical moments. In particular, the expected value �(t) ≡ 〈�r(t)〉 (14) provides the best estimate of the spreading of the heterogeneous fronts and is adopted as the definition of mixing length. We refer to Ref. [22] for a detailed discussion. 5 Scaling Laws and Mixing Regimes Both nonlinearity of the flow equations and permeability heterogeneity can cause dispersive mixing of the fluid transport. To understand their combined effect on the mixing process, we first clarify the effect of each separately. 5.1 The effect of nonlinearity A linear stability analysis of the water/oil front in a homogeneous porous medium is carried out in, e.g., Ref. [6]. This analysis reveals the existence of two distinct flow regimes, characterized in terms of the frontal mobility ratio Λ = λ(s−)/λ(s+): The front is stable when Λ < 1 and unstable when Λ > 1. For our choice of constitutive functions (f and λ), Λ < 1 (respectively Λ > 1) when μo/μw < 2.657 (resp. μo/μw > 2.657). The linearized stability theory predicts that for Λ > 1 small front perturbations will grow exponentially. The long-time nonlinear development of front perturbations is studied numerically in Ref. [21]. In brief, to a first approximation, in an unstable flow process �(t) = O(t), as t → ∞. (15) 5.2 The effect of heterogeneity We consider linear displacement (tracer flow) in heterogeneous porous media. In tracer flows, the two fluids are miscible and have equal viscosities. In this case, the frontal 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 63 mobility ratio Λ ≡ 1 and the flow is neutrally stable. As a result, the mixing in tracer flows is solely the effect of velocity dispersion caused by heterogeneity. For tracer flow in fractal permeability fields, the theory of fluid mixing discussed in Ref. [26] provides the following scaling laws: �(t) = O(tγ ), with γ = max { 1 2 , 1 2 + 1 + β 2 } . (16) Two qualitatively distinct regimes are seen. If β ≤ −1, then γ = 1/2, and the mixing process is Fickian. If −1 < β < 0, then γ > 1/2, and the mixing is anomalous (i.e., the diffusivity increases with time or, equivalently, with travel distance). These scaling laws were consistently derived from leading order results of primitive and renormalized perturbation theory. They were also verified for large perturbation parameters by numerical simulations in Ref. [22]. 5.3 The combined effect of nonlinearity and heterogeneity The analysis of fluid mixing dynamics in heterogeneous, two-phase flow is via compu- tational experiments. The computations involved (1) the construction of the various ensembles of random permeability fields, characterized by different values of the ex- ponent β and the coefficient of variation CVk, and (2) the subsequent determination of the resultant mixing region for the fluid flows through these ensembles. Both nu- merical and statistical convergence were verified by successive mesh refinements and increase in ensemble size, respectively. All results reported are for a flow region with Lx/Ly = 4. All log-permeability fields were drawn from a Gaussian distribution with the power law covariance (7). Two permeability exponents, β = −∞, −0.5, and sev- eral values for the permeability coefficient of variation, CVk, were considered. We use the latter as a dimensionless measure of heterogeneity. The rationale for the choice of exponents is the desire to explore the effect of short length (β = −∞) and long length scale (β = −0.5) heterogeneities in two-phase flow dispersion. Table 1 compiles our results. This table displays the values of M , CVk, and β used in each study, and the corresponding mixing regimes. The regimes are classified according to the asymptotic scaling of �(t) = O(tδ), as t → ∞: N U (Nonlinear Unstable) if δ = 1; N S (Nonlinear Stable) if δ = 0; L (Linear) if δ = γ, the scaling exponent for the linear flow problem given in (16); L+ (Superlinear) if γ < δ < 1; L− (Sublinear) if 0 < δ < γ. Table 1: Mixing regimes for two-phase flow CVk/M 0.5 1 2.657 5 10 20 β 0.54 N S L− L N U N U N U −∞ 1.33 L− L L+ N U N U −∞ 2.93 L− L L L+ N U −∞ 0.49 L− L− L+ N U N U −0.5 0.99 L− L N U N U N U −0.5 1.83 L− L L N U −0.5 64 Frederico Furtado and Felipe Pereira 6, 4(2004) Our findings can be summarized as follows: • Distinct mixing regimes occur depending on whether one of the driving mecha- nisms (nonlinearity and heterogeneity) dominates the mixing dynamics or not. • If nonlinearity dominates, then two distinct mixing regimes are possible, accord- ing to whether the nonlinear effects are stabilizing (N S regime) or destabilizing (N U regime). • If heterogeneity dominates, then the linear regime L is observed. In this case, mixing regions in linear and nonlinear flows grow at essentially identical asymp- totic rates. • Intermediate mixing regimes (sublinear and superlinear) occur when nonlinear- ity and heterogeneity compete for the dominance of the mixing dynamics. Figure 4 is a pictorial representation of these results in terms of a “phase-diagram” in the M vs. CVk plane. The results suggest that for a given value of β there is a critical value for M (Mcritical in the figure). At nearby values, the two-phase flow system behaves according to the several mixing regimes discussed above. The location of Mcritical seems to depend on β, being smaller for larger values of β (c.f. Table 1). Figure 4: A pictorial representation of the distinct mixing regimes that result from the inter- play between nonlinearity and heterogeneity: Nonlinearity controlled regimes N S (nonlinear stable) and N U (nonlinear unstable); heterogeneity controlled regime L (linear); intermedi- ate regimes L− (sublinear) and L+ (superlinear). Figure 5 illustrates the crossover from heterogeneity controlled to nonlinearity con- trolled mixing as the heterogeneity strength is varied while the nonlinearity strength is kept fixed. The figure refers to a study in which β = −∞, M = 5 were kept fixed and the strength of heterogeneity was increased by increasing CVk. 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 65 Figure 5: Large time mixing length curves for two-phase flows (M = 5.0) in β = −∞ permeability fields corresponding to distinct mixing regimes. From top to bottom the CVk values are 0.54, 1.33 and 2.93 and the mixing regimes are nonlinear unstable N U , superlinear L+ and linear L, respectively. 66 Frederico Furtado and Felipe Pereira 6, 4(2004) Each picture in this figure shows the log-log plot of the mixing length �(t) as a function of the travel distance (the irregular curve). Also plotted are two straight lines with vertical positions fixed to agree with the last plotted point on the mixing length curves. The slopes of these lines are 0.5, the asymptotic scaling exponent for the linear mixing regime L for this value of β, and 1, the asymptotic scaling exponent of the nonlinear unstable mixing regime N U . It is clear that as the CVk is increased the mixing length curves bend continuously from the N U regime towards the L regime. We close this section with a few remarks on the scale up problem. For weak heterogeneities the mixing process is ‘convective’, with the mixing region growing lin- early as a function of time. This behavior is typical of mixing processes driven by the nonlinear instability of the flow and should be properly modeled with a hyperbolic renormalization of the flow equations. However, for sufficiently strong heterogeneities with short correlation length, the mixing process is Fickian and a dispersive renormal- ization, like the one developed recently for linear flow [26], should perhaps be more appropriate. 6 Computational Solutions: The Decoupling Hy- pothesis The studies [30, 31, 45, 32, 14] assume that the transport and saturation equations are not coupled. Under this assumption distinct theoretical developments for the scale up of two-phase flow systems have been reported in these papers. Here we test such assumption by computing directly the mixing length growth in a fixed (in time) velocity field. We consider the nonlinear unstable N U mixing regime characterized by weak heterogeneities. New ensemble averages have been performed using the sets of parameters (CVk,M ,β) displayed in Table 1 which produced the nonlinear unstable mixing regime N U . In these new ensemble average studies the velocity field was computed at the first time step of the simulations (with the given intial condition for the saturation) and was kept fixed in time. Thus, the computations were restricted to the solution of the saturation equation (2). All the new ensemble average studies fall in the linear mixing regime L. The asymptotic scaling exponents for the linear mixing regimes which we found were 0.5 for the ensemble averages with β = −∞ and 0.75 for the ensemble averages having β = −0.5 for the Hurst exponent. Figures 6 and 7 illustrate our new results for β = −∞ and β = −0.5, respectively. Consider Figure 6 first. It refers to the case CVk = 0.54 and M = 5.0. This figure illustrates the change in the mixing regime as the velocity field is kept fixed in time. Each picture in this figure shows the log-log plot of the mixing length �(t) as a function of the travel distance (the irregular curve). Also plotted are two straight lines with vertical positions fixed to agree with the last plotted point on the mixing length curves. The slopes of these lines are 0.5, the asymptotic scaling exponent for the linear mixing regime L for β = −∞, and 1, the asymptotic scaling exponent of the 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 67 Figure 6: The effect of decoupling saturation and pressure equations. If the coupling is taken into account we have determined that the mixing regime is the nonlinear unstable N U (curve labeled with variable velocity) for the parameters CVk = 0.54, M = 5.0 and β = −∞. Once the equations are decoupled, the computed mixing regime is the linear L (curve labeled with constant velocity). Clearly the dark solid line would lead to incorrect predictions for the behavior of the mixing process as a function of time. Figure 7: The effect of decoupling saturation and pressure equations. If the coupling is taken into account we have determined that the mixing regime is the nonlinear unstable N U (curve labelled with variable velocity) for the parameters CVk = 0.99, M = 10.0 and β = −0.5. Once the equations are decoupled, the computed mixing regime is the linear L (curve labelled with constant velocity). Predictions made with the constant (in time) velocity field would be very inaccurate. 68 Frederico Furtado and Felipe Pereira 6, 4(2004) nonlinear unstable mixing regime N U . It is clear that as the velocity field is kept fixed in time the mixing length curves change from the N U regime to the L regime. It is also clear from these results that the new mixing length curves would lead to incorrect predictions for the behavior of the mixing process as a function of time. Figure 7 refers to a study for the parameters CVk = 0.99 and M = 10.0. The same behavior for the mixing length growth described by Figure 6 was observed here. As the velocity field is kept fixed the nonlinear unstable mixing regime N U changes to the linear mixing regime L. Consequently predictions made with the constant (in time) velocity field would be very inaccurate. The authors are now investigating the validity of the decoupling hypothesis in the description the linear mixing regime L. 7 Conclusions In this paper we reviewed our recent effort to characterize fluid mixing regimes in multiphase flow in multiscale heterogeneous porous formations. We believe such a characterization is an essential step in the scientific program leading to an improved understanding of the scale up problem for porous media flow. Our recent results on the characterization of mixing regimes for two-phase, immis- cible, incompressible displacement in petroleum reservoirs have been used to investi- gate a hypothesis that has been used in the literature in the development of theoretical predictions for the scale up problem. According to these studies the saturation and pressure equations are not coupled. Highly resolved numerical simulations were conducted to study this hypothesis and we determined that for the mixing regime characterized by weak heterogeneities (or, equivalently, by strong nonlinearities in the governing equations) it would lead to incorrect predictions. The decoupling hypothesis for the linear mixing regime, characterized by strong heterogeneities, is currently being investigated. The effect of gravity and the presence of stochastic layers in the geology are im- portant problems to be considered within the computational stochastic approach de- veloped by the authors. References [1] Y. Amirat, K. Hamdache, and A. Ziani, Homogénéisation d’équations hyper- boliques du premier ordre et application aux écoulements miscibles en milieu poreux, Ann. Inst. H. Poincaré 6 (1989) 397. [2] J. A. Aronovitz and D. R. Nelson, Anomalous diffusion in steady fluid flow through a porous medium, Phys. Rev. A 30 (1984) 1948–1954. [3] L. An, J. Glimm, D. H. Sharp, and Q. Zhang. Scale up of flow in porous media. In A. P. Bourgeat, C. Carasso, S. Luckhaus, and A. Mikelić, editors, Mathematical 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 69 Modelling of Flow Through Porous Media, pages 26–44. World Scientific, New Jersey, 1995. [4] A. Arya, T. Hewett, R. Larson, L. Lake, Dispersion and reservoir heterogeneity, SPE Res. Eng. J. 3 (1988) 139-148. [5] J. W. Barker and S. Thibeau. A critical review of the use of pseudorelative per- meabilities for upscaling. SPE Reservoir Engineering, 138–143, May 1997. [6] A. J. Chorin. The instability of fronts in a porous medium. Comm. Math. Phys., 91 (1983) 103–116. [7] M. A. Christie. Upscaling for reservoir simulations. JPT J. Pet. Technol., 48 (1996) 1004-1010. [8] J. Cushman. Physics of Fluids in Hierarchical Porous Media: Angstrons to Miles, Cambridge, in press, 1997. [9] V. Cvetkovic and G. Dagan. Reactive transport and immiscible flow in geological media. II. Applications. Proc. R. Soc. Lond., A 452 (1996) 303–328. [10] G. Dagan. Dispersion of a passive solute in non-ergodic transport by steady ve- locity fields in heterogeneous formations. J. Fluid Mech. 233 (1991) 1281-1290. [11] G. Dagan. The significance of heterogeneity of evolving scales to transport in porous formations. Water Resour. Res. 30 (1994) 3327-3336. [12] G. Dagan, Theory of solute transport by groundwater, Annu. Rev. Fluid Mech. 19 (1987) 183–215. [13] G. Dagan, Flow and Transport in Porous Formations, Springer-Verlag, New York, 1989. [14] G. Dagan and V. Cvetkovic, Reactive transport and immiscible flow in geological media. I. General Theory. Proc. R. Soc. Lond., A 452 (1996) 285–301. [15] J. Douglas, Jr., F. Furtado, and F. Pereira. On the numerical simulation of waterflooding of heterogeneous petroleum reservoirs. Computational Geosciences 1(2) (1997) 155-190. [16] L. J. Durlofsky. Use of higher moments for the description of upscaled, process independent relative permeabilities. SPE 37987, 1997. [17] Y. 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. [18] F. Furtado, J. Glimm, B. Lindquist, and F. Pereira. Multi-length scale calcula- tions of mixing length growth in tracer floods. In F. Kovarik, editor, Proceedings of the Emerging Technologies Conference, pages 251–259. Institute for Improved Oil Recovery, U. Houston, Houston, TX, 1990. 70 Frederico Furtado and Felipe Pereira 6, 4(2004) [19] F. Furtado, J. Glimm, B. Lindquist, and F. Pereira, Characterization of mixing length growth for flow in heterogeneous porous media, SPE 21233 (1991). [20] F. Furtado, J. Glimm, B. Lindquist, F. Pereira, and Q. Zhang, Time dependent anomalous diffusion for flow in multi-fractal porous media, Proceedings, Work- shop on Numerical Methods for the Simulation of Multiphase and Complex Flow, ed. by T.M.M. Verheggan, Lecture Notes in Physics, vol. 398, Springer Verlag, New York, 1992, 211–220. [21] F. Furtado and F. Pereira. Scaling analysis for two-phase, immiscible flow in heterogeneous media. Computational and Applied Mathematics, 17 (1998) 233- 262. [22] F. Furtado and F. Pereira. Crossover from nonlinearity controlled to heterogeneity controlled mixing in two-phase porous media flows. Computational Geosciences 2003. To appear. [23] L. W. Gelhar, Stochastic subsurface hydrology from theory to applications, Water Resour. Res. 22 (1986) 135S–145S. [24] J. Glimm, H. Kim, D. Sharp, and T. Wallstrom. A stochastic analysis of the scale up problem for flow in porous media. Computational and Applied Mathematics, 17 (1998) 67–79. [25] J. Glimm, B. Lindquist, F. Pereira, and R. Peierls, The fractal hypothesis and anomalous diffusion, Computational and Applied Mathematics, 11 (1992) 189– 207. [26] J. Glimm, B. Lindquist, F. Pereira, and Q. Zhang. A theory of macrodispersion for the scale up problem. Transport in Porous Media, 13 (1993) 97–122. [27] J. Glimm and D. H. Sharp, A random field model for anomalous diffusion in heterogeneous porous media, J. Stat. Phys., 62 (1991) 415–424. [28] J. Glimm and D.H. Sharp. Stochastic partial differential equations: selected ap- plications in continuum physics. In R. A. Carmona and B. L. Rozovskii, editors, Stochastic Partial Differential Equations: Six Perspectives, Mathematical Sur- veys and Monographs, American Mathematical Society, Providence, 1997. [29] T. A. Hewett, Fractal distributions of reservoir heterogeneity and their influence on fluid transport, SPE 15386, 1986. [30] K. D. Jarman, Jr. Stochastic immiscible flow with moment equations. Ph.D. The- sis, University of Colorado, Boulder, 2000. [31] K. D. Jarman and T. F. Russell. Analysis of 1-D moment equations for immiscible flow. Contemporary Mathematics. To appear. [32] P. Langlo and M. Espedal. Macrodispersion for two-phase, immiscible flow in porous media. Adv. Water Resour., 17 (1994) 297–316. 6, 4(2004) On the Scale Up Problem for Two-Phase Flow in Petroleum Reservoirs 71 [33] R. Lenormand. Determining flow equations from stochastic properties of a per- meability field: the MHD model. SPE Journal, 179–190, June 1996. [34] H. Nessyahu and E. Tadmor. Non-oscillatory central differencing for hyperbolic conservation laws. J. Comp. Phys. 87(2) (1990) 408-463. [35] S. P. Neuman and Y. K. Zhang, A quasi-linear theory of non-Fickian subsurface dispersion, Water Resour. Res. 26 (1990) 887–902. [36] S. Neuman, Universal Scaling of Hydraulic Conductivities and Dispersivities in Geologic Media, Water Resour. Res. 26 (1990) 1794-1758. [37] D. W. Peaceman, Fundamentals of Numerical Reservoir Simulation, Elsevier, New York, 1977. [38] F. Pereira, Stochastic geology and porous media flow: theory and simulations, State University of New York at Stony Brook, Ph.D. Thesis, 1992. [39] I. Gel’fand and N. Vilenkin. Generalized Functions , IV (English translation), Academic Press, New York, 1964. [40] T. Wallstrom, S. Hou, M. A. Christie, L. J. Durlofsky, and D. H. Sharp. Accurate scale up of two-phase flow using renormalization and nonuniform coarsening. Computational Geoscience 3 (1999) 69-87. [41] T. Wallstrom, S. Hou, M. A. Christie, L. J. Durlofsky, and D. H. Sharp. Applica- tion of a new two-phase upscaling technique to realistic reservoir cross sections. Proceedings of the SPE 15th Symposium on Reservoir Simulation, pages 451-462, 1999. SPE 51939. [42] T. Wallstrom, S. Hou, M. A. Christie, L. J. Durlofsky, D. H. Sharp, and Q. Zhou. Effective medium boundary conditions for upscaling relative permeabilities. Transport in Porous Media, 2000. Accepted for publication. [43] J. R. Waggoner, J. L. Castillo, and L. W. Lake. Simulation of EOR processes in stochastically generated permeable media. SPE Formation Evaluation, June 1992, 173-180. SPE 21237. [44] C. L. Winter, C. M. Neuman, and S. P. Newman, A perturbation expansion for diffusion in a random velocity field, SIAM J. Appl. Math. 44 (1984) 411–424. [45] D. Zhang, L. Li, and H. A. Tchelepi. Stochastic formulation for uncertainty analysis of two-phase flow in heterogeneous reservoirs. SPE Journal 5 (2000) 60–70. [46] Q. Zhang. A multi-length scale theory of the anomalous mixing length growth for tracer flow in heterogeneous porous media. J. Stat. Phys., 66 (1991) 485–501. [47] Q. Zhang, The asymptotic scaling behavior of mixing induced by a random velocity field, Adv. Appl. Math., 16 (1995) 23–58. 72 Frederico Furtado and Felipe Pereira 6, 4(2004) [48] Q. Zhang, The transient behavior of mixing induced by a random velocity field, Water Resour. Res. 31, (1995) 577–591. [49] J. Zhu. A numerical study of stochastic dispersion for two-phase immiscible flow in porous media. M.Sc. Thesis, University of Wyoming, 2001.