Journal articles on the topic 'Multiphase Upscaling'

Summary Representing the complete spectrum of fine-scale stratigraphic details in full-field dynamic models of geologically complex clastic reservoirs is beyond the limits of existing computational capabilities. A quasiglobal multiphase upscaling method—the regional-scale multiphase upscaling (RMU) method—is developed, in which the dynamic effects of subgrid-scale (typically subseismic) nonlocal stratigraphic reservoir elements (e.g., channels, lobes, sand bars, and shale drapes) are captured by means of pseudofunctions for flow simulation. Unlike conventional dynamic multiphase upscaling methods, the RMU method does not require fine-resolution reservoir-scale simulations. Rather, it relies on intermediate-scale sector-model simulations for pseudoization. The intermediate scale, also referred to as the regional scale, is defined as the spatial scale at which the global multiphase flow effects of nonlocal stratigraphic elements can be approximated by fine-resolution flow simulations with reasonable accuracy. During the pseudoization process, dynamic multiphase flow responses of coarse regional-scale sector models are calibrated against those stemming from their corresponding fine-resolution parent models. Each regional-scale sector model is simulated only once at the fine geologic resolution. The process involves automatic determination and subsequent modification of the parameters that describe rock relative permeability and capillary pressure functions. Coarse regional-scale models are simulated a few times until a reasonable match between their coarse- and fine-resolution dynamic responses can be attained. The parameter-estimation step of the pseudoization process is performed by use of a very efficient constrained nonlinear optimization algorithm. The RMU method is evaluated in two proof-of-concept numerical examples involving a plethora of turbidite stratigraphic architectures. The method yields simulation results that are always more accurate than conventionally upscaled coarse-resolution model predictions. Incorporating geologically based pseudofunctions into otherwise simple coarse-resolution full-field reservoir models reduces the simulation cycle time significantly and improves the accuracy of production forecasts. The RMU method typically delivers two to three orders of magnitude in speed up of flow simulations.

Summary Carbonate fractured reservoirs introduce a tremendous challenge to the upscaling of both single- and multiphase flow. The complexity comes from both heterogeneous matrix and fracture systems in which the separation of scales is very difficult. The mathematical upscaling techniques, derived from representative elementary volume (REV), must therefore be replaced by a more realistic geology-based approach. In the case of multiphase flow, an evaluation of the main forces acting during oil recovery must also be performed. A matrix-sector model from a highly heterogeneous carbonate reservoir is linked to different fracture realizations in dual-continuum simulations. An integrated iterative workflow between the geology-based static modeling and the dynamic simulations is used to investigate the effect of fracture heterogeneity on multiphase fluid flow. Heterogeneities at various scales (i.e., diffuse fractures and subseismic faults) are considered. The diffuse-fracture model is built on the basis of facies and porosity from the matrix model together with core data, image-log data, and data from outcrop-analogs. Because of poor seismic data, the subseismic-fault model is mainly conceptual and is based on the analysis of outcrop-analog data. Fluid-flow simulations are run for both single-phase and multiphase flow and gas and water injections. A better understanding of fractured-reservoirs behavior is achieved by incorporating realistic fracture heterogeneity into the geological model and analyzing the dynamic impact of fractures at various scales. In the case of diffuse fractures, the heterogeneity effect can be captured in the upscaled model. The subseismic faults, however, must be explicitly represented, unless the sigma (shape) factor is included in the upscaling process. A local grid-refinement approach is applied to demonstrate explicit fractures in large-scale simulation grids. This study provides guidelines on how to effectively scale up a heterogeneous fracture model and still capture the heterogeneous flow behavior.

Summary Although numerous upscaling techniques are reported in the literature, efficiently computing reasonably accurate equivalent rock properties from geological data at fine scale remains difficult. This is especially true for facies with multiple lithologies under multiphase-flow conditions. Because of the nature of multiscale heterogeneity inherent in petroleum reservoirs, the equivalent rock and flow properties will vary with the scales of heterogeneity. Therefore, upscaled properties under multiphase-flow conditions cannot be estimated without reference to the absolute scales of heterogeneity. Wavelet analysis is a multiresolution framework; thus, it is well suited for upscaling rock and flow properties in a multiscale heterogeneous reservoir. The compact support property of the wavelet transform assures efficient computation. Choice of regularity provides a flexible way to control the smoothness of the resulting upscaling properties. In this study, we developed a new procedure to improve significantly the computational efficiency and accuracy of upscaling for generating equivalent rock and rock-fluid properties under various geological and flow conditions based on multiresolution analysis of wavelet transforms. Additionally, we explored a wavelet reconstruction method to provide a basis for downsampling fine-scale rock property fields from information at various levels of coarser scale. The beauty of the method is that because the equivalent properties at different length scales are computed recursively, the interdependent influences of the heterogeneities on the scales are included effectively. We demonstrate the method by successfully applying it to upscale interbedset and interfacies reservoir properties of Almond formation outcrops under multiphase-flow conditions. Introduction It is well known that petroleum reservoirs are inherently heterogeneous, flow performance of reservoirs is controlled by variability in reservoir properties at various scales, and the dominant scale effects will vary with the production processes involved. Limited by the computational power or based on the consideration of the dominant scale on the production scheme involved, the flow simulation is usually carried on a scale that is much coarser than the scale of the direct or indirect measurements. Therefore, one needs first to reconstruct adequately the fine-scale variability mapping of the reservoir properties from sparsely distributed, multiscale samples, and then properly coarsen, or upscale the fine-scale variability. Roughly speaking, upscaling is considered a procedure to transfer the rock and rock-fluid properties carried in the finer scale to appropriate coarser scale. Upscaling reservoir properties is complicated and challenging, especially for multiple lithologies and multiphase-flow conditions. Unlike the case of single-phase flow, where only absolute permeability affects the steady-state flow, where only absolute permeability affects the steady-state flow behavior, the flow properties for each phase and phase pressure differ within each geological facies, and, probably more important, the rock-fluid property contrasts between the facies also play a vital role for multiple lithologies and multiphase-flow situations. Therefore, besides effective absolute permeability, as in the case of single-phase flow, one needs to determine two additional equivalent properties: effective relative permeabilities and effective capillary pressure for the coarse-scale grid model. Such equivalent properties are flow-regime sensitive, making the problem more complicated. It is well known that flow regime in a multiphase flow is determined by the relative ratio among three major driving forces: viscous force, gravitational force, and capillary pressure. Such a relative ratio among various driving forces will depend on flow velocity, scale of heterogeneity, size of gridblock, and fluid properties. Considering that flow velocity and heterogeneity in a reservoir are spatially and/or temporally varying the equivalent relative permeability and capillary pressure generally will be spatially and temporally dependent and cannot be constructed without reference to geological structure and dominant flow regime. Although the literature reports numerous upscaling techniques,1–5 efficiently computing reasonably accurate equivalent rock properties from the data at finer scale remains challenging. We lack a method that can be used consistently for both single and multiphase-flow conditions. With multiresolution analysis properties, compact support, orthogonality, and localized transforms in both space and frequency domains, wavelet analysis provides a consistent, efficient, and accurate way for upscaling reservoir properties. In this study, we developed a decomposition and reconstruction procedure for rock property fields in light of Mallat's work.6 Application of the developed procedure to absolute permeability fields presents promising results. When coarsening fine-grid permeability to more accessible scales, although the number of grids are reduced substantially, the main trends or characters of the permeability fields are reasonably well preserved. The flow behavior from the original fine-grid permeability field is reproduced perfectly by the upscaled permeability field. The developed method also can reconstruct the fine-scale properties adequately from coarser-scale information. To check the consistency of the proposed procedure, the same framework was extended further to upscale rock-fluid properties under multiple lithologies and multiphase-flow conditions in addition to rock properties. To verify the applicability of the newly developed upscaling procedure within a realistic geological context, we constructed testing problems with measurements from the Almond formation. Theoretical Background Wavelet Transforms. Wavelet transforms, like the extensively used Fourier transform, are kinds of linear integral transforms that are represented as a convolution of a given function, f(x), with a kernel function,Equation 1 Here, ? is a scale parameter, x is a location parameter, and the family of functions, ? ?x(u), are called wavelet functions. There are two major differences between the Fourier transform and Wavelet transforms. First of all, instead of one parameter transforms, as Eq. 1 shows. In this way, Wavelet transforms give not only information about the content of a function in the frequency domain, but also information about the location of these frequencies in the spatial domain. Wavelet Transforms. Wavelet transforms, like the extensively used Fourier transform, are kinds of linear integral transforms that are represented as a convolution of a given function, f(x), with a kernel function,Equation 1 Here, ? is a scale parameter, x is a location parameter, and the family of functions, ? ?x(u), are called wavelet functions. There are two major differences between the Fourier transform and Wavelet transforms. First of all, instead of one parameter transforms, as Eq. 1 shows. In this way, Wavelet transforms give not only information about the content of a function in the frequency domain, but also information about the location of these frequencies in the spatial domain.

Summary Geologists often generate highly heterogeneous descriptions of reservoirs, containing complex structures which are likely to give rise to very tortuous flow paths. However, these models contain too many grid cells for multiphase flow simulation, and the number of cells must be reduced by upscaling for reservoir simulation. Conventional upscaling methods often have difficulty in the representation of tortuous flow paths, mainly because of the inappropriate assumptions concerning the boundary conditions. An accurate and practical upscaling method is therefore required to preserve the flow features caused by highly heterogeneous fine scale geological description. In this paper, the problems encountered in routinely used upscaling approaches are outlined, and a more accurate and practical way of performing upscaling is proposed. The new upscaling method, Well Drive Upscaling (WDU), employs the wells and the actual reservoir boundary conditions (e.g., faults and physical boundaries of the geological model). The main advantage of this method is that the dominant flow paths can be preserved, and thus the geological knowledge can be assimilated appropriately. The new method has firstly been applied to a synthetic model with a tortuous channel, and is shown to have significant improvement over the traditional approach. The sensitivity study on the scale-up factor using a benchmark model shows the advantage of the method with various scale-up factors. The method was then applied to a model of a field in the central North Sea, which involves three-phase flow. In the cases studied, the WDU method produced a comparable result to the dynamic Pore Volume Weighted approach, which involves running the fine grid simulation and computing appropriate relative permeabilities and interblock transmissibilities. The new method makes the upscaling process practical, and our tests show it to be more accurate than traditional methods. Introduction The heterogeneity observed in a field is generally high and the geological structures therein can be complex. From a geological point of view, it would be ideal to represent each facies boundary, both vertically and horizontally, by a gridblock boundary (Mallet 1997; Deutsch and Tran 2002). Also, if distinct layering exists within a genetic unit, a further split into subunits is also desirable. In practice, reservoir models are usually created at the scale of meters or less vertically and 100 meters or less areally [and each block itself may have involved small-scale upscaling (Pickup et al. 2005)]. In many cases, detailed reservoir modeling for a highly heterogeneous reservoir may result in a large number of grid cells (e.g., 106 grid cells or more). This large number of grid cells prohibits direct simulation of the reservoir, especially for a very heterogeneous reservoir model. This is because, apart from the limitation of computational power, the high level of heterogeneity often makes it difficult to obtain a converged solution. The problem becomes more severe when simulations involve three-phase flow. In order to perform reservoir simulation on a highly heterogeneous geological model within a reasonable time frame, we have to apply appropriate upscaling techniques to reduce the number of grid cells so as to speed up the reservoir simulation and thus field development planning process. Although a number of upscaling methods have been developed in the past a few decades (Pickup et al. 2005; Christie 1996, 2001), they are often not satisfactory and have been discussed in a number of critical reviews (Barker and Thibeau 1997; Farmer 2002). The main conflict in the application of the current upscaling techniques lies in the balance of the accuracy and practicality of the methods. There are two main problems that cause the conflict. The first is the problem of using inappropriate boundary conditions in single-phase upscaling, which is likely to reduce the accuracy, and the second problem is the impracticality of the dynamic two-phase upscaling which should (in theory) be more accurate. Details of the these methods have been discussed in a number of reviews on upscaling (Christie 1996, 2001; Barker and Thibeau 1997; Farmer 2002; Renard and de Marsily 1997), so a complete review of upscaling methods will not be presented here. However, we outline one of the commonly used methods: the pressure solution method for upscaling single-phase flow. In this method, a single-phase pressure solve is carried out in each coarse cell in turn, and Darcy's law is used to calculate the effective permeability tensor (Christie 1996). In order to solve the pressure equation, boundary conditions must be applied to each cell. (This is referred to as the local upscaling method.) A typical example is the no-flow, or constant pressure boundary condition, where the pressure is fixed at either end of the region of a coarse block, and no flow is allowed through the sides. Other boundary conditions include linear pressure and periodic boundary conditions (Farmer 2002). Such boundary conditions, however, may differ significantly from the actual boundary conditions within a heterogeneous fine-scale model (Chen et al. 2003; Zhang 2006). A highly heterogeneous reservoir model often produces tortuous flow paths, and it is difficult to generalize a flow pattern on the boundaries of a coarse block and apply to all the coarse blocks in a reservoir model. The flow paths for a coarse model may be completely different from the original fine-scale geological model response when inappropriate boundary conditions are applied (i.e., the effect of geological structure may be lost). A multiphase flow simulation from such a coarse model will not honor the small-scale geological structure either, even if we ignore the error caused by multiphase flow effects in the upscaling process.

Summary Upscaling is often applied to coarsen detailed geological reservoir descriptions to sizes that can be accommodated by flow simulators. Adaptive local-global upscaling is a new and accurate methodology that incorporates global coarse-scale flow information into the boundary conditions used to compute upscaled quantities (e.g., coarse-scale transmissibilities). The procedure is iterated until a self-consistent solution is obtained. In this work, we extend this approach to 3D systems and introduce and evaluate procedures to decrease the computational demands of the method. This includes the use of purely local upscaling calculations for the initial estimation of coarse-scale transmissibilities and the use of reduced border regions during the iterations. This is shown to decrease the computational requirements of the reduced procedure significantly relative to the full methodology, while impacting the accuracy very little. The performance of the adaptive local-global upscaling technique is evaluated for three different heterogeneous reservoir descriptions. The method is shown to provide a high degree of accuracy for relevant flow quantities. In addition, it is shown to be less computationally demanding and significantly more accurate than some existing extended local upscaling procedures. Introduction Fine-scale heterogeneity can have a significant impact on reservoir performance. Because it is usually not feasible to simulate directly on the detailed geocellular model, some type of upscaling is often applied to generate the simulation model from the geological description. Here, we focus on the upscaling of single-phase flow parameters, particularly absolute permeability. The algorithms we consider can provide either coarse-scale permeability, designated k*, or coarse-scale transmissibility, designated T* . It is important to emphasize that the accurate upscaling of permeability (which can be studied within the context of single-phase flow) is essential for the development of accurate coarse models of two-phase or multiphase flow. Thus the applicability of the methods developed here is very broad and includes all types of displacement processes.

Fractured reservoirs are distributed widely over the world, and describing fluid flow in fractures is an important and challenging topic in research. Discrete fracture modeling (DFM) and equivalent continuum modeling are two principal methods used to model fluid flow through fractured rocks. In this paper, a novel method, embedded discrete fracture modeling (EDFM), is developed to compute equivalent permeability in fractured reservoirs. This paper begins with an introduction on EDFM. Then, the paper describes an upscaling procedure to calculate equivalent permeability. Following this, the paper carries out a series of simulations to compare the computation cost between DFM and EDFM. In addition, the method is verified by embedded discrete fracture modeling and fine grid methods, and grid-block and multiphase flow are studied to prove the feasibility of the method. Finally, the upscaling procedure is applied to a three-dimensional case in order to study performance for a gas injection problem. This study is the first to use embedded discrete fracture modeling to compute equivalent permeability for fractured reservoirs. This paper also provides a detailed comparison and discussion on embedded discrete fracture modeling and discrete fracture modeling in the context of equivalent permeability computation with a single-phase model. Most importantly, this study addresses whether this novel method can be used in multiphase flow in a reservoir with fractures.

Summary Enhanced-oil-recovery (EOR) processes involve complex flow, transport, and thermodynamic interactions; as a result, compositional simulation is necessary for accurate representation of the physics. Flow simulation of compositional systems with high-resolution reservoir models is computationally intensive because of the large number of unknowns and the strong nonlinear interactions. Thus, there is a great need for upscaling methods of compositional processes. The complex multiscale interactions between the phase behavior and the heterogeneities lie at the core of the difficulty in constructing consistent upscaling procedures. We use a mass-conservative formulation and introduce upscaled phase-molar-mobility functions for coarse-scale modeling of multiphase flow. These upscaled flow functions account for the subgrid effects caused by the absolute permeability and relative permeability variations, as well as the effects of compressibility. Upscaling of the phase behavior is performed as follows. We assume that instantaneous thermodynamic equilibrium is valid on the fine scale, and we derive coarse-scale equations in which the phase behavior may not necessarily be at equilibrium. The upscaled thermodynamic functions, which represent differences in the component fugacities, are used to account for the nonequilibrium effects on the coarse scale. We demonstrate that the upscaled phase-behavior functions transform the equilibrium phase space on the fine scale to a region of similar shape, but with tilted tie-lines on the coarse space. The numerical framework uses K-values that depend on the orientation of the tie-lines in the new nonequilibrium phase space and the sign of upscaled thermodynamic functions. The proposed methodology is applied to challenging gas-injection problems with large numbers of components and highly heterogeneous permeability fields. The K-value-based coarse-scale operator produces results that are in good agreement with the fine-scale solutions for the quantities of interest, including the component overall compositions and saturation distributions.

Summary One difficulty in fracture upscaling for field-scale dual-porosity reservoir simulation is the determination of equivalent gridblock fracture permeability, which depends on the type of boundary conditions imposed on the discrete-fracture-network (DFN) simulation. Actually, classical upscaling procedures usually are based on linearly varying pressure boundary conditions, which cannot capture the near-well flow behavior. As a result, the well productivity calculated by a dual-porosity flow simulator can be very different from that calculated on a DFN model. This paper proposes a near-well fracture-upscaling procedure based on the geological DFN model to improve the accuracy of well productivity in fractured-reservoir simulators. This procedure enables us to represent the actual flow through the fractures and the exchanges between matrix and fractures in the well vicinity. On the basis of the computed near-well flow pattern, equivalent fracture transmissibilities as well as numerical well indices are determined and assigned to the gridblocks of the dual-porosity reservoir simulator. The reliability and necessity of using the near-well upscaling procedure are demonstrated by examples. Introduction Advanced characterization methodologies are now able to provide realistic models of geological fracture networks (Cacas et al. 2001). In addition, production logging and transient well tests can be simulated with DFN models to validate the geological fracture-network geometry and calibrate the hydraulic properties of fractures (Sarda et al. 2002). However, because of computational limitations, the complex geological DFN model cannot be used straightforwardly to simulate a multiphase-flow production scenario at field scale (Bourbiaux et al. 2002). For such simulations, a dual-porosity reservoir simulator is typically used. The dual-porosity reservoir model, using large gridblocks to discretize the whole reservoir, is a conceptual representation of the actual geology of the fractured medium. The flow properties of the fracture network are then homogenized on gridblocks through upscaling procedures. The upscaling of fracture properties is the problem of translating the geological and hydraulic description of fracture networks into reservoir-simulation parameters. The dual-porosity model requires the determination of equivalent fracture permeability and equivalent matrix-block dimensions or shape factors (Bourbiaux et al. 1997; Sarda et al. 1997). This paper discusses methodologies for upscaling the permeability of a fracture network, especially in the vicinity of the well. Upscaling of fracture permeability has been studied extensively. The commonly used method is numerical, based on flow simulation on a model of the actual fracture network with specific boundary conditions to compute an equivalent gridblock permeability (Sarda et al. 1997). Other methods were also developed; for example, Oda (1985) proposed an analytical equation to calculate the fracture-permeability tensor, and Lough et al. (1997) presented an approach using the boundary-element method, which integrates the contribution of matrix in the equivalent permeability of the fractured medium. When using a numerical approach to determine the equivalent permeability of a fracture network, the upscaled result depends on the type of boundary conditions imposed in the flow simulation. Actually, classical upscaling procedures are usually based on flow simulation in a parallelepipedic model with linear-type pressure boundary conditions, which cannot capture the near-well flow behavior. As a result, the well productivity calculated by a dual-porosity flow simulator can be very different from that calculated on a near-wellbore DFN model.

Summary Upscaling is often needed in reservoir simulation to coarsen highly detailed geological descriptions. Most existing upscaling procedures aim to reproduce fine-scale results for a particular geological model (realization). In this work, we develop and test a new approach, ensemble-level upscaling, for efficiently generating upscaled two-phase flow parameters (e.g., upscaled relative permeabilities) for multiple geological realizations. The ensemble-level upscaling approach aims to achieve agreement between the fine- and coarse-scale flow models at the ensemble level, rather than realization-by-realization agreement, as is the intent of existing upscaling techniques. For this purpose, flow-based upscaling calculations are combined with a statistical procedure based on a cluster analysis. This approach allows us to compute numerically the upscaled two-phase flow functions for only a small fraction of the coarse blocks. For the majority of blocks, these functions are estimated statistically on the basis of single-phase velocity information (attributes), determined when the upscaled single-phase parameters are calculated. The procedure is designed to maintain close correspondence between the cumulative distribution functions (CDFs) for the numerically computed and statistically estimated two-phase flow functions. We apply the method to 2D synthetic models of multiple realizations for uncertainty quantification. Models with different geological heterogeneity and fluid-mobility ratios are considered. It is shown that the method consistently corrects the biases evident in primitive coarse-scale predictions and can capture the ensemble statistics (e.g., P50, P10, P90) of the fine-scale results almost as accurately as the full flow-based upscaling procedures but with much less computational effort. The overall approach is flexible and can be used with any combination of upscaling procedures. Introduction In recent years, a wide variety of upscaling procedures has been developed and applied. These techniques generally take as their starting point a fine-scale geological model of the subsurface. The intent is then to generate a coarser model, which retains the geological realism of the underlying fine-scale description, for use in flow simulation. Though model sizes can vary substantially depending on the application, typical fine-scale geocellular models may contain 107 to 108 cells, while typical simulation models may contain 104 to 106 blocks. Recent reviews and assessments (e.g., Barker and Thibeau 1997; Barker and Dupouy 1999; Farmer 2002; Darman et al. 2002; Gerritsen and Durlofsky 2005; Chen 2005) describe and apply a variety of upscaling techniques. These procedures can be categorized in different ways. One important distinction is in terms of the coarse-scale parameters that are computed by a particular method. Specifically, a technique that generates only upscaled single-phase parameters (permeability or transmissibility) can be classified as a single-phase upscaling procedure even though it may be applied to two- or three-phase flow problems. A method that additionally generates upscaled relative permeability functions is termed a two-phase upscaling procedure. Another way to distinguish upscaling procedures is according to the problem solved to determine the coarse-scale parameters. In particular, methods may be classified as local, extended local, quasiglobal, or global in order of increasing computational effort, depending on the problem solved in the upscaling computations. In general, two-phase upscaling methods are more computationally expensive than single-phase upscaling procedures, as a time-dependent two-phase flow problem must be solved in this case. The appropriate upscaling procedure for any particular problem depends on the required level of accuracy and the degree of coarsening. For example, for permeability fields characterized by two-point geostatistics (variogram-based models), with only a moderate degree of coarsening, the use of local single-phase upscaling procedures, possibly coupled with nonuniform gridding, may provide acceptable coarse models. For more challenging cases, however, such as channelized systems characterized by multipoint geostatistics and high degrees of upscaling, extended local or (quasi) global single-phase upscaling coupled with two-phase upscaling may be necessary. In recent work (Chen and Durlofsky 2006b), we introduced an upscaling procedure that combines quasiglobal single-phase upscaling, which was accomplished through a local-global procedure, with a specialized two-phase upscaling. The technique was shown to provide reasonable degrees of accuracy for challenging problems, though it was observed that the speedups between fine-grid simulation and the upscaling plus coarse-scale simulations were not that dramatic (e.g., approximately a factor of 4 to 10). Speedups will be much more substantial if the model is simulated many times, because the computation time required for the two-phase upscaling calculations is large compared to the coarse-grid simulations. It would, however, still be useful to accelerate these upscaling computations. This is particularly desirable in cases with substantial uncertainty in the underlying geological model, in which case many realizations (or scenarios) are to be simulated. In such cases, realization-by-realization agreement between fine and coarse models is less essential. Rather, what is required in this case is agreement of a statistical nature, such as agreement in the CDFs (e.g., the P10, P50, P90 predictions) for relevant production quantities such as cumulative oil recovered or net present value. The required level of accuracy of the upscaling, on the realization-by-realization basis, could be slightly less for such cases, though the method should be unbiased. The intent of this paper is to develop and test procedures for substantially accelerating two-phase upscaling procedures for cases in which many realizations are to be considered. Toward this goal, we couple upscaling with statistical estimation techniques. Several statistical techniques were considered, though the best performance was achieved using K-means clustering. Application of this approach allows us to compute upscaled two-phase functions through full-flow simulation for only a small fraction of the coarse-scale blocks. For the rest of the blocks, these functions are estimated statistically on the basis of velocity information (attributes) computed during the single-phase upscaling. The overall method can be used with any combination of single-phase and two-phase upscaling procedures and is shown to provide a high level of accuracy in the statistical sense described above for example cases involving different heterogeneity models. There has been very little research reported on the development of upscaling procedures for multiple permeability realizations. Previous researchers considered related problems involving the handling of upscaled multiphase flow parameters (e.g., the grouping of pseudorelative permeabilities). Dupouy et al. (1998) applied a statistical procedure to group the numerically computed global pseudorelative permeabilities to reduce the number of pseudofunctions used in flow simulation. Their work did not involve the estimation of upscaled relative permeabilities, though they noted that such an approach would be useful in practice because it would reduce the number of pseudofunctions to be numerically computed. Christie and Clifford (1998) suggested an a priori approach to grouping upscaled parameters for compositional simulation. They used the concept of tracer-breakthrough curves to represent coarse-scale blocks, and applied K-means clustering analysis to group the upscaled functions. Neither of these studies, however, considered upscaling over multiple reservoir models and the associated assessment of uncertainty for fine-scale predictions. Our work here is also related to previous studies on error modeling of coarse-scale simulation models (Omre and Lødøen 2004, Lødøen et al. 2004), though the approaches are quite different. In the error modeling studies, some fine-scale calibration runs were required to model upscaling error and correct the bias in the coarse-scale simulation results, while our approach here estimates the upscaled flow parameters directly. The statistical estimation procedure (based on cluster analysis) used here can be viewed as a proxy or surrogate method that avoids the need to numerically generate upscaled two-phase parameters. In this sense, any proxy can be applied in the procedure. Statistical clustering approaches are used in many applications and have been applied recently in reservoir engineering as proxies for simulations in genetic algorithm-based optimization (Artus et al. 2006). The outline of this paper is as follows: We first provide the governing equations and a brief overview of the relevant upscaling procedures. Next, we describe and illustrate the ensemble-level upscaling approach based on clustering to estimate statistically the upscaled two-phase flow functions. This is followed by extensive numerical results for a variety of 2D systems. We conclude with a discussion and summary.

Summary Simulation of an in-situ combustion (ISC) process was performed for a fractured system at core and matrix-block scales. The aim of this work was: (1) To predict the ISC extinction/propagation condition(s), (2) understand the mechanism of oil recovery, and (3) provide some guidelines for ISC upscaling for a fractured system. The study was based on a fine-grid, single-porosity, multiphase, and multicomponent simulation using a thermal reservoir simulator. First, the simulator was validated for 1D combustion using the corresponding analytical solutions. 2D combustion was validated using experimental results available in the literature. It was found that the grid size should not be larger than the combustion-zone thickness in order for the results to be independent of grid size. ISC in the fractured system was strongly dependent on the oxygen diffusion coefficient, while the matrix permeability played an important role in oil production. The effect of each production mechanism was studied separately whenever it was possible. Oil production is governed mainly by oil drainage because of gravity force, which is enhanced by viscosity reduction; possible pressure-gradient generation in the ISC process seems to have a minor effect. The nature (oil-production rate, saturations distribution, shape of the combustion front) of ISC at core scale was different from that in a single block with surrounding fracture. The important characteristics of different zones (i.e., combustion, coke, and oil zones) at block scale were studied, and some preliminary guidelines for upscaling are presented.

We use the method of density functional hydrodynamics (DFH) to model compositional multiphase flows in natural cores at the pore-scale. In previous publications the authors demonstrated that DFH covers many diverse pore-scale phenomena, starting from those inherent in RCA and SCAL measurements, and extending to much more complex EOR processes. We perform the pore-scale modelling of multiphase flow scenarios by means of the direct hydrodynamic (DHD) simulator, which is a numerical implementation of the DFH. In the present work, we consider the problem of numerical modelling of fluid transport in pore systems with voids and channels when the range of pore sizes exceed several orders of magnitude. Such situations are well known for carbonate reservoirs, where narrow pore channels of micrometer range can coexist and interconnect with vugs of millimeter or centimeter range. In such multiscale systems one cannot use the standard DFH approach for pore-scale modeling, primarily because the needed increase in scanning resolution that is required to resolve small pores adequately, leads to a field of view reduction that compromises the representation of large pores. In order to address this challenge, we suggest a novel approach, in which transport in small-size pores is described by an upscaled effective model, while the transport in large pores is still described by the DFH. The upscaled effective model is derived from the exact DFH equations using asymptotic expansion in respect to small-size characterization parameter. This effective model retains the properties of DFH like chemical and multiphase transport, thus making it applicable to the same range of phenomena as DFH itself. The model is based on the concept that the transport is driven by gradients of chemical potentials of the components present in the mixture. This is a significant generalization of the Darcy transport model since the proposed new model incorporates diffusion transport in addition to the usual pressure-driven transport. In the present work we provide several multiphase transport numerical examples including: a) upscaling to chemical potential drive (CPD) model, b) combined modeling of large pores by DFH and small pores by CPD.

Water infiltration and unsaturated flow through heterogeneous soil control the distribution of soil moisture in the vadose zone and the dynamics of groundwater recharge, providing the link between climate, biogeochemical soil processes and vegetation dynamics. Infiltration into dry soil is hydrodynamically unstable, leading to preferential flow through narrow wet regions (fingers). In this paper we use numerical simulation to study the interplay between fingering instabilities and soil heterogeneity during water infiltration. We consider soil with heterogeneous intrinsic permeability. Permeabilities are random, with point Gaussian statistics, and vary smoothly in space due to spatial correlation. The key research question is whether the presence of moderate or strong heterogeneity overwhelms the fingering instability, recovering the simple stable displacement patterns predicted by most simplified model of infiltration currently used in hydrological models from the Darcy to the basin scales. We perform detailed simulations of constant-rate infiltration into soils with isotropic and anisotropic intrinsic permeability fields. Our results demonstrate that soil heterogeneity does not suppress fingering instabilities, but it rather enhances its effect of preferential flow and channeling. Fingering patterns strongly depend on soil structure, in particular the correlation length and anisotropy of the permeability field. While the finger size and flow dynamics are only slightly controlled by correlation length in isotropic fields, layering leads to significant finger meandering and bulging, changing arrival times and wetting efficiencies. Fingering and soil heterogeneity need to be considered when upscaling the constitutive relationships of multiphase flow in porous media (relative permeability and water retention curve) from the finger to field and basin scales. While relative permeabilities remain unchanged upon upscaling for stable displacements, the inefficient wetting due to fingering leads to relative permeabilities at the field scale that are significantly different from those at the Darcy scale. These effective relative permeability functions also depend, although less strongly, on heterogeneity and soil structure.

This paper (SPE 52637) was revised for publication from paper SPE 38002, first presented at the 1997 SPE Reservoir Simulation Symposium, Dallas, 8-11 June. Original manuscript received for review 1 July 1997. Revised manuscript received 5 August 1998. Paper peer approved 3 September 1998. Summary The gridblock permeabilities used in reservoir simulation are commonly determined through the upscaling of a fine scale geostatistical reservoir description. Though it is well established that permeabilities computed in this manner are, in general, full tensor quantities, most finite difference reservoir simulators still treat permeability as a diagonal tensor. In this paper, we implement a capability to handle full tensor permeabilities in a general purpose finite difference simulator and apply this capability to the modeling of several complex geological systems. We formulate a flux continuous approach for the pressure equation by use of a method analogous to that of previous researchers (Edwards and Rogers; Aavatsmark et al.), consider methods for upwinding in multiphase flow problems, and additionally discuss some relevant implementation and reservoir characterization issues. The accuracy of the finite difference formulation, assessed through comparisons to an accurate finite element approach, is shown to be generally good, particularly for immiscible displacements in heterogeneous systems. The formulation is then applied to the simulation of upscaled descriptions of several geologically complex reservoirs involving crossbedding and extensive fracturing. The method performs quite well for these systems and is shown to capture the effects of the underlying geology accurately. Finally, the significant errors that can be incurred through inaccurate representation of the full permeability tensor are demonstrated for several cases. P. 567

Summary The control-volume discrete-fracture (CVDF) model is extended to incorporate heterogeneity in rock and in rock-fluid properties. A novel algorithm is proposed to model strong water-wetting with zero capillary pressure in the fractures. The extended method is used to simulate:oil production in a layered faulted reservoir,laboratory displacement tests in a stack of matrix blocks with a large contrast in fracture and matrix capillary pressure functions, andwater injection in 2D and 3D fractured media with mixed-wettability state. Our results show that the algorithm is suitable for the simulation of water injection in heterogeneous porous media both in water-wet and mixed-wettability states. The novel approach with zero fracture capillary and nonzero matrix capillary pressure allows the proper prediction of sharp fronts in the fractures. Introduction This work is focused on the numerical treatment of two main physical aspects of multiphase flow in fractured porous media: heterogeneity in rock-fluid properties and reservoir wettability. In a previous work (Monteagudo and Firoozabadi 2004), a CVDF method was used to discretize the system of equations governing water injection in fractured media with strong-water-wettability state and homogeneous matrix and rock-fluid properties. The method was restricted to a finite contrast in matrix-fracture capillary pressure. In this work, we extend the CVDF model for simulation of water injection in fractured media comprised of heterogeneous rocks and wettability conditions from strong-water-wetting to mixed-wetting conditions. We also present a formulation for infinite contrast in capillary pressures of matrix and fractures (zero capillary pressure in the fracture and finite capillary pressure in the matrix). The control volume (CV) method, first proposed by Baliga and Patankar (1980), is a finite-volume formulation over dual cells (CVs) of a Delaunay mesh. It is locally conservative and suited for unstructured grids. It has been widely employed for the simulation of multiphase flow in porous media (Monteagudo and Firoozabadi 2004; Verma 1996; Helmig 1997; Helmig and Huber 1998; Bastian et al. 2000; Geiger et al. 2003) and the convergence of the method for two-phase immiscible flow in porous medium has already been proved (Michel 2003). Numerical treatment of heterogeneity in the framework of the CV method has been extensively studied in the past (Edwards 2002; Edwards and Rogers 1998; Prevost 2000; Aavatsmark et al. 1998a, b). Nevertheless, those works have focused on absolute permeability heterogeneity and anisotropy in single-phase flow. The main concern in those works is the use of full tensor permeability and the accurate generation of streamlines (required by the streamline numerical method). It is well known that the standard CV method produces inaccurate velocity fields around the interfaces of heterogeneous media as the contrast in permeability is increased (Durlofsky 1994). In the standard CV method, Delaunay triangles are locally homogeneous and the polygonal CV cell may be heterogeneous (see Fig. 1a). For accurate streamlines, several authors (Verma 1996; Edwards 2002; Edwards and Rogers 1998; Prevost 2000; Aavatsmark et al. 1998a) have proposed that the polygonal CV cell must be locally homogeneous, implying heterogeneous Delaunay triangles (see Fig. 1b). The latter configuration, however, generates additional problems in the simulation of multiphase flow in porous media. Basically, from mesh generation standpoint, it may not be possible to generate an unstructured mesh where the boundaries of the CV median-dual cell conform to heterogeneous interfaces in the domain. Conforming mesh is important for the discrete-fracture approach. Therefore, it would be necessary to first generate a standard CV cell mesh, and later a homogenization procedure would be required to obtain CV cells with constant permeability. The homogenization or upscaling of permeability is somehow possible, but the same is not true for rock-fluid properties; most challenging is capillary pressure with different endpoints. Therefore, the approach with the homogeneous CV cell may be suitable for single-phase simulation where rock-fluid interactions are not part of the problem. However, rock-fluid interactions have to be taken into account for simulation of multiphase flow in fractured porous medium. Frequently, capillary pressure is disregarded in two-phase flow simulations; however, capillary pressure is of importance for simulation of multiphase flow in fractured porous media (Monteagudo and Firoozabadi 2004; Karimi-Fard and Firoozabadi 2003). Predictions of flow pattern and oil recovery may be severely affected if capillary pressure effect is neglected.

Abstract. There is a great interest and demand for green-type energy storage in Sweden for both short- and long-term (hours, days, weeks and seasons) periods. While there are a number of approaches proposed (e.g., compressed air, geothermal and thermal), only a few have commercially been demonstrated through upscaling projects. Among these, the thermal energy storage (TES) that stores energy (excess heat or cold) in fluids is particularly interesting. The excess energy can be stored underground in excavated caverns and used for large district heating and cooling purposes as well as for balancing and regulating electrical energy in power grids. For an upscaling underground TES project within the Tornquist suture zone of Scania in the southwest of Sweden, three high-resolution seismic profiles, each approximately 1 km long, were acquired. Geologically, the site sits within the southern margin of the Romeleåsen fault zone in the Sorgenfrei–Tornquist Zone (STZ), where dolerite dyke swarms of Carboniferous–Permian age are observed striking in the SE–NW direction for hundreds of kilometers both on land and in offshore seismic and magnetic data (from Scania to Midland Valley in the UK). These dykes, 10–50 m thick, in the nearby quarries (within both Precambrian gneiss and quartzite) express themselves mostly in a subvertical manner. They can therefore act as a good water/fluid barrier, which can be an important geological factor for any TES site. For the data acquisition, combined cabled and wireless recorders were used to provide continuity on both sides of a major road running in the middle of the study area. Bedrock depressions are clearly depicted in the tomograms, suggesting the possibility of zones of weaknesses, highly fractured and/or weathered, in the bedrock and confirmed in several places by follow-up boreholes. Several steeply dipping (60–65°) reflections were imaged down to 400 m depth and interpreted to originate from dolerite dykes. This interpretation is based on their orientations, strong amplitudes, regular occurrences and correlation with downhole logging data. In addition, groundwater flow measurements within the unconsolidated sediments and in bedrock suggest steeply dipping structures are the dominant factor in directing water mainly along a SE–NW trend, which is consistent with the strike of the dyke swarm within the STZ. To provide further insight on the origin of the reflections, even the historical crustal-scale offshore BABEL (Baltic and Bothnian Echoes from the Lithosphere) lines (A-AA-AB) were revisited. Clear multiphase faults and signs of intrusions or melt source in the lower crust are observed, as well as a Moho step across the Tornquist zone. Overall, we favor that the reflections are of dolerite origin and their dip component (i.e., not subvertical) may imply a Precambrian basement (and dykes) tilting, block rotation, towards the NE as a result of the Romeleåsen reverse faulting. In terms of thermal storage, these dykes then may be encountered during the excavation of the site and can complicate underground water flow should they be used as a fluid barrier in case of leakage.

Summary Advancements in horizontal-well drilling and multistage hydraulic fracturing have enabled economically viable gas production from tight formations. Reservoir-simulation models play an important role in the production forecasting and field-development planning. To enhance their predictive capabilities and to capture the uncertainties in model parameters, one should calibrate stochastic reservoir models to both geologic and flow observations. In this paper, a novel approach to characterization and history matching of hydrocarbon production from a hydraulic-fractured shale is presented. This new methodology includes generating multiple discrete-fracture-network (DFN) models, upscaling the models for numerical multiphase-flow simulation, and updating the DFN-model parameters with dynamic-flow responses. First, measurements from hydraulic-fracture treatment, petrophysical interpretation, and in-situ stress data are used to estimate the initial probability distribution of hydraulic-fracture and induced-microfracture parameters, and multiple initial DFN models are generated. Next, the DFN models are upscaled into an equivalent continuum dual-porosity model with analytical techniques. The upscaled models are subjected to the flow simulation, and their production performances are compared with the actual responses. Finally, an assisted-history-matching algorithm is implemented to assess the uncertainties of the DFN-model parameters. Hydraulic-fracture parameters including half-length and transmissivity are updated, and the length, transmissivity, intensity, and spatial distribution of the induced fractures are also estimated. The proposed methodology is applied to facilitate characterization of fracture parameters of a multifractured shale-gas well in the Horn River basin. Fracture parameters and stimulated reservoir volume (SRV) derived from the updated DFN models are in agreement with estimates from microseismic interpretation and rate-transient analysis. The key advantage of this integrated assisted-history-matching approach is that uncertainties in fracture parameters are represented by the multiple equally probable DFN models and their upscaled flow-simulation models, which honor the hard data and match the dynamic production history. This work highlights the significance of uncertainties in SRV and hydraulic-fracture parameters. It also provides insight into the value of microseismic data when integrated into a rigorous production-history-matching work flow.

Summary We use upscaling through homogenization to predict oil recovery from fractured reservoirs consisting of matrix columns, also called vertically fractured reservoirs (VFRs), for a variety of conditions. The upscaled VFR model overcomes limitations of the dual-porosity model, including the use of a shape factor. The purpose of this paper is to investigate three main physical aspects of multiphase flow in fractured reservoirs: reservoir wettability, viscosity ratio, and heterogeneity in rock/fluid properties. The main characteristic that determines reservoir behavior is the Péclet number that expresses the ratio of the average imbibition time divided by the residence time of the fluids in the fractures. The second characteristic dimensionless number is the gravity number. Upscaled VFR simulations, aimed at studying the mentioned features, add new insights. First, we discuss the results at low Péclet numbers. For only small gravity numbers, the effect of contact angle, delay time for the nonequilibrium capillary effect, the heterogeneity of the matrix-column size, and the matrix permeability can be ignored without appreciable loss of accuracy. The ultimate oil recovery for mixed-wet VFRs is approximately equal to the Amott index, and the oil production does not depend on the absolute value of the phase viscosity but on viscosity ratio. However, large gravity numbers enhance underriding, aggravated by large contact angles, longer delay times, and higher viscosity ratios. Layering can lead to an improvement or deterioration, depending on the fracture aperture and permeability distribution. At low Péclet numbers, the fractured reservoir behaves very similarly to a conventional reservoir and depends largely on the viscosity ratio and the gravity number. At high Péclet numbers, after water breakthrough, the oil recovery appears to be proportional to the cosine of the contact angle and inversely proportional to the sum of the oil and water viscosity. In addition, the mixed-wetting effect is more pronounced; there are significant influences of delay time (nonequilibrium effects), matrix permeability, matrix-column size, and the column-size distribution on oil recovery. At low gravity numbers and an effective length/thickness ratio larger than 10, the oil recovery is independent of the vertical-fracture-aperture distribution. For the same amount of injected water, the recovery at low Péclet numbers is larger than the recovery at high Péclet numbers.

Summary Multiscale methods have been developed for accurate and efficient numerical solution of flow problems in large-scale heterogeneous reservoirs. A scalable and extendible Operator-Based Multiscale Method (OBMM) is described here. OBMM is cast as a general algebraic framework. It is natural and convenient to incorporate more physics in OBMM for multiscale computation. In OBMM, two operators are constructed: prolongation and restriction. The prolongation operator is constructed by assembling the multiscale basis functions. The specific form of the restriction operator depends on the coarse-scale discretization formulation (e.g., finitevolume or finite-element). The coarse-scale pressure equation is obtained algebraically by applying the prolongation and restriction operators to the fine-scale flow equations. Solving the coarse-scale equation results in a high-quality coarse-scale pressure. The finescale pressure can be reconstructed by applying the prolongation operator to the coarse-scale pressure. A conservative fine-scale velocity field is then reconstructed to solve the transport (saturation) equation. We describe the OBMM approach for multiscale modeling of compressible multiphase flow. We show that extension from incompressible to compressible flows is straightforward. No special treatment for compressibility is required. The efficiency of multiscale formulations over standard fine-scale methods is retained by OBMM. The accuracy of OBMM is demonstrated using several numerical examples including a challenging depletion problem in a strongly heterogeneous permeability field (SPE 10). Introduction The accuracy of simulating subsurface flow relies strongly on the detailed geologic description of the porous formation. Formation properties such as porosity and permeability typically vary over many scales. As a result, it is not unusual for a detailed geologic description to require 107-108 grid cells. However, this level of resolution is far beyond the computational capability of state-of-the-art reservoir simulators (106 grid cells). Moreover, in many applications, large numbers of reservoir simulations are performed (e.g., history matching, sensitivity analysis and stochastic simulation). Thus, it is necessary to have an efficient and accurate computational method to study these highly detailed models. Multiscale formulations are very promising due to their ability to resolve fine-scale information accurately without direct solution of the global fine-scale equations. Recently, there has been increasing interest in multiscale methods. Hou and Wu (1997) proposed a multiscale finite-element method (MsFEM) that captures the fine-scale information by constructing special basis functions within each element. However, the reconstructed fine-scale velocity is not conservative. Later, Chen and Hou (2003) proposed a conservative mixed finite-element multiscale method. Another multiscale mixed finite element method was presented by Arbogast (2002) and Arbogast and Bryant (2002). Numerical Green functions were used to resolve the fine-scale information, which are then coupled with coarse-scale operators to obtain the global solution. Aarnes (2004) proposed a modified mixed finite-element method, which constructs special basis functions sensitive to the nature of the elliptic problem. Chen et al. (2003) developed a local-global upscaling method by extracting local boundary conditions from a global solution, and then constructing coarse-scale system from local solutions. All these methods considered incompressible flow in heterogeneous porous media where the pressure equation is elliptic. A multiscale finite-volume method (MsFVM) was proposed by Jenny et al. (2003, 2004, 2006) for heterogeneous elliptic problems. They employed two sets of basis functions--dual and primal. The dual basis functions are identical to those of Hou and Wu (1997), while the primal basis functions are obtained by solving local elliptic problems with Neumann boundary conditions calculated from the dual basis functions. Existing multiscale methods (Aarnes 2004; Arbogast 2002; Chen and Hou 2003; Hou and Wu 1997; Jenny et al. 2003) deal with the incompressible flow problem only. However, compressibility will be significant if a gas phase is present. Gas has a large compressibility, which is a strong function of pressure. Therefore, there can be significant spatial compressibility variations in the reservoir, and this is a challenge for multiscale modeling. Very recently, Lunati and Jenny (2006) considered compressible multiphase flow in the framework of MsFVM. They proposed three models to account for the effects of compressibility. Using those models, compressibility effects were represented in the coarse-scale equations and the reconstructed fine-scale fluxes according to the magnitude of compressibility. Motivated to construct a flexible algebraic multiscale framework that can deal with compressible multiphase flow in highly detailed heterogeneous models, we developed an operator-based multiscale method (OBMM). The OBMM algorithm is composed of four steps:constructing the prolongation and restriction operators,assembling and solving the coarse-scale pressure equations,reconstructing the fine-scale pressure and velocity fields, andsolving the fine-scale transport equations. OBMM is a general algebraic multiscale framework for compressible multiphase flow. This algebraic framework can also be extended naturally from structured to unstructured grid. Moreover, the OBMM approach may be used to employ multiscale solution strategies in existing simulators with a relatively small investment.

Summary Reservoir heterogeneity can be detrimental to the success of surfactant/polymer enhanced-oil-recovery (EOR) processes. Therefore, it is important to evaluate the effect of uncertainty in reservoir heterogeneity on the performance of surfactant/polymer EOR. Usually, a Monte Carlo sampling approach is used, in which a number of stochastic reservoir-model realizations are generated and then numerical simulation is performed to obtain a certain objective function, such as the recovery factor. However, Monte Carlo simulation (MCS) has a slow convergence rate and requires a large number of samples to produce accurate results. This can be computationally expensive when using large complex reservoir models. This study applies a multiscale approach to improve the efficiency of uncertainty quantification. This method is known as the multilevel Monte Carlo (MLMC) method. This method comprises performing a small number of expensive simulations on the fine-scale model and a large number of less-expensive simulations on coarser upscaled models, and then combining the results to produce the quantities of interest. The purpose of this method is to reduce computational cost while maintaining the accuracy of the fine-scale model. The results of this approach are compared with a reference MCS, assuming a large number of simulations on the fine-scale model. Other advantages of the MLMC method are its nonintrusiveness and its scalability to incorporate an increasing number of uncertainties. This study uses the MLMC method to efficiently quantify the effect of uncertainty in heterogeneity on the recovery factor of a chemical EOR process, specifically surfactant/polymer flooding. The permeability field is assumed to be the random input. This method is first demonstrated by use of a Gaussian 3D reservoir model. Different coarsening algorithms are used and compared, such as the renormalization method and the pressure-solver method (PSM). The results are compared with running Monte Carlo for the fine-scale model while equating the computational cost for the MLMC method. Both of these results are then compared with the reference case, which uses a large number of runs of the fine-scale model. The method is then extended to a channelized non-Gaussian generated 3D reservoir model incorporating multiphase upscaling The results show that it is possible to robustly quantify spatial uncertainty for a surfactant/polymer EOR process while greatly reducing the computational requirement, up to two orders of magnitude compared with traditional Monte Carlo for both the Gaussian and non-Gaussian reservoir models. The method can be easily extended to other EOR processes to quantify spatial uncertainty, such as carbon dioxide (CO2) EOR. Other possible extensions of this method are also discussed.

Summary The modern focused ion beam-scanning electron microscopy (FIB-SEM) allows imaging of nanoporous tight reservoir-rock samples in 3D at a resolution up to 3 nm/voxel. Correct porosity determination from FIB-SEM images requires fast and robust segmentation. However, the quality and efficient segmentation of FIB-SEM images is still a complicated and challenging task. Typically, a trained operator spends days or weeks in subjective and semimanual labeling of a single FIB-SEM data set. The presence of FIB-SEM artifacts, such as porebacks, requires developing a new methodology for efficient image segmentation. We have developed a method for simplification of multimodal segmentation of FIB-SEM data sets using machine-learning (ML)-based techniques. We study a collection of rock samples formed according to the petrophysical interpretation of well logs from a complex tight gas reservoir rock of the Berezov Formation (West Siberia, Russia). The core samples were passed through a multiscale imaging workflow for pore-space-structure upscaling from nanometer to log scale. FIB-SEM imaging resolved the finest scale using a dual-beam analytical system. Image segmentation used an architecture derived from a convolutional neural network (CNN) in the DeepUNet (Ronneberger et al. 2015) configuration. We implemented the solution in the Pytorch® (Facebook, Inc., Menlo Park, California, USA) framework in a Linux environment. Computation exploited a high-performance computing system. The acquired data included three 3D FIB-SEM data sets with a physical size of approximately 20 × 15 × 25 µm with a voxel size of 5 nm. A professional geologist manually segmented (labeled) a fraction of slices. We split the labeled slices into training, validation, and test data. We then augmented the training data to increase its size. The developed CNN delivered promising results. The model performed automatic segmentation with the following average quality indicators according to test data: accuracy of 86.66%, precision of 54.93%, recall of 83.76%, and F1 score of 55.10%. We achieved a significant boost in segmentation speed of 14.5 megapixel (MP)/min. Compared with 0.18 to 1.45 MP/min for manual labeling, this yielded an efficiency increase of at least 10 times. The presented research work improves the quality of quantitative petrophysical characterization of complex reservoir rocks using digital rock imaging. The development allows the multiphase segmentation of 3D FIB-SEM data complicated with artifacts. It delivers correct and precise pore-space segmentation, resulting in little turn-around-time saving and increased porosity-data quality. Although image segmentation using CNNs is mainstream in the modern ML world, it is an emerging novel approach for reservoir-characterization tasks.

Summary Modern streamline-based reservoir simulators are able to account for actual field conditions such as 3D multiphase flow effects, reservoir heterogeneity, gravity, and changing well conditions. A streamline simulator was used to model four field cases, with approximately 400 wells and 150,000 gridblocks. History-match run times were approximately 1 CPU hour per run, with the final history matches completed in approximately 1 month per field. In all field cases, a high percentage of wells were history matched within the first two to three runs. Streamline simulation not only enables a rapid turnaround time for studies, but it also serves as a different tool in resolving each of the studied fields' unique characteristics. The primary reasons for faster history matching of permeability fields using 3D streamline technology as compared to conventional finite-difference (FD) techniques are as follows: Streamlines clearly identify which producer-injector pairs communicate strongly (flow visualization). Streamlines allow the use of a very large number of wells, thereby substantially reducing the uncertainty associated with outer-boundary conditions. Streamline flow paths indicate that idealized drainage patterns do not exist in real fields. It is therefore unrealistic to extract symmetric elements out of a full field. The speed and efficiency of the method allows the solution of fine-scale and/or full-field models with hundreds of wells. The streamline simulator honors the historical total fluid injection and production volumes exactly because there are no drawdown constraints for incompressible problems. The technology allows for easy identification of regions that require modifications to achieve a history match. Streamlines provide new flow information (i.e., well connectivity, drainage volumes, and well allocation factors) that cannot be derived from conventional simulation methods. Introduction In the past, streamline-based flow simulation was quite limited in its application to field data. Emanuel and Milliken1 showed how hybrid streamtube models were used to history match field data rapidly to arrive at both an updated geologic model and a current oil-saturation distribution for input to FD simulations. FD simulators were then used in forecast mode. Recent advances in streamline-based flow simulators have overcome many of the limitations of previous streamline and streamtube methods.2-6 Streamline-based simulators are now fully 3D and account for multiphase gravity and fluid mobility effects as well as compressibility effects. Another key improvement is that the simulator can now account for changing well conditions due to rate changes, infill drilling, producer-injector conversions, and well abandonments. With advances in streamline methods, the technique is rapidly becoming a common tool to assist in the modeling and forecasting of field cases. As this technology has matured, it is becoming available to a larger group of engineers and is no longer confined to research centers. Published case studies using streamline simulators are now appearing from a broad distribution of sources.7–12 Because of the increasing interest in this technology, our first intent in this paper is to outline a methodology for where and how streamline-based simulation fits in the reservoir engineering toolbox. Our second objective is to provide insight into why we think the method works so well in some cases. Finally, we will demonstrate the application of the technology to everyday field situations useful to mainstream exploitation or reservoir engineers, as opposed to specialized or research applications. The Streamline Simulation Method For a more detailed mathematical description of the streamline method, please refer to the Appendix and subsequent references. In brief, the streamline simulation method solves a 3D problem by decoupling it into a series of 1D problems, each one solved along a streamline. Unlike FD simulation, streamline simulation relies on transporting fluids along a dynamically changing streamline- based flow grid, as opposed to the underlying Cartesian grid. The result is that large timestep sizes can be taken without numerical instabilities, giving the streamline method a near-linear scaling in terms of CPU efficiency vs. model size.6 For very large models, streamline-based simulators can be one to two orders of magnitude faster than FD methods. The timestep size in streamline methods is not limited by a classic grid throughput (CFL) condition but by how far fluids can be transported along the current streamline grid before the streamlines need to be updated. Factors that influence this limit include nonlinear effects like mobility, gravity, and well rate changes.5 In real field displacements, historical well effects have a far greater impact on streamline-pattern changes than do mobility and gravity. Thus, the key is determining how much historical data can be upscaled without significantly impacting simulation results. For all cases considered here, 1-year timestep sizes were more than adequate to capture changes in historical data, gravity, and mobility effects. It is worth noting that upscaling historical data also would benefit run times for FD simulations. Where possible, both SL and FD methods would then require similar simulation times. However, only for very coarse grids and specific problems is it possible to take 1-year timestep sizes with FD methods. As the grid becomes finer, CFL limitations begin to dictate the timestep size, which is much smaller than is necessary to honor nonlinearities. This is why streamline methods exhibit larger speed-up factors over FD methods as the number of grid cells increases.

Summary The shape factor concept, originally introduced by Barenblatt in 1960, provides an elegant and powerful upscaling method for fractured reservoir simulation. The shape factor determines the fluid and heat transfer between matrix and fractures when there is a difference in pressure or temperature between matrix blocks and the surrounding fractures. An appropriate specification of the shape factor is therefore critical for accurate modeling. Since its introduction, many different values for the shape factor have been proposed in the literature, among which the well-known Warren-Root and Kazemi shape factors. The aim of this paper is to show that the selection of the appropriate shape factor should not only depend on the "shape" and dimensions of matrix blocks, but should also take into consideration the character of the dominant underlying physical recovery mechanisms. We will show that by taking into account the dominant physical recovery mechanism, the apparent discrepancies in the shape factor values reported in the literature can be overcome. We derive a general expression for the shape factor that not only captures existing shape factor expressions, but also allows extensions to recovery mechanisms requiring a dual permeability approach. The paper is organized as follows. First, we briefly review the shape factors presented in the literature. We then derive the general expression for the (single-phase) matrix-fracture shape factor. Subsequently, we analytically derive a new shape factor that captures the transient in pressure/temperature diffusion processes. To compare and contrast the impact of the various shape factors, we consider three cases of increasing complexity. First, we consider pressure/temperature diffusion in a single 1D matrix block following a step change in the boundary conditions. Next, we consider isothermal gas/oil gravity drainage from a homogeneous stack. We compare fine-grid single-porosity simulations (in which the matrix is finely gridded and in which the fractures are explicitly represented) with coarse-grid dual-permeability simulations (in which the matrix-fracture interaction is modeled by shape factors). In the third step, we consider gas-oil gravity drainage of the same stack model, but now under steam injection. In this case, steam is injected at the top, and oil recovered from the base of the fracture system. Again, we compare fine-grid single-porosity simulations with coarse-grid dual-permeability simulations. We show that in this case, the constant (asymptotic) shape factor provides a good approximation to the heating of the stack. We will show, however, that with a constant (time-independent) shape factor, the initial fast heating of the matrix blocks cannot be captured. We show that the new transient shape factor, however, enables coarse-grid dual-permeability modeling of thermal recovery processes such that they reproduce fine-grid results. Introduction The modeling of matrix-fracture interaction using shape factors has been an active area of research for over 40 years now, and has attracted considerable attention both in the context of single- and multi-phase matrix-fracture modeling (Barenblatt et al. 1960; Warren and Root 1963; Kazemi et al. 1976; Thomas et al. 1983; Coats 1989; Ueda et al. 1989; Zimmerman et al. 1993a; Chang 1993; Lim and Aziz 1995; Gilman and Kazemi 1983; Beckner et al. 1987, 1988; Rossen and Shen 1989; Bech et al. 1991; Bourbiaux et al. 1999). In their 1960 landmark paper, Barenblatt et al. introduced the shape factor concept to model the (single-phase) fluid transfer between matrix and fractures (1960). The central idea of Barenblatt et al. was not to study the behavior of individual matrix blocks and their surrounding fractures, but instead to introduce two abstract interacting media: one medium, the "matrix," in which the physical matrix blocks are lumped, and one medium, the "fractures," in which the fractures are lumped. Whenever a pressure difference exists between the matrix and the fractures, a fluid flow between the media will occur. The shape factor is then defined by the following relation, which ties the (single-phase) matrix-fracture fluid flow to the instantaneous pressure difference between matrix and fractures:q = s (km / µ) V (p*m - pf), ....[ EQ. 1 ] where V denotes the volume of the matrix block. In 1963, Warren and Root used Barenblatt's shape factor concept in the context of well-testing using dual porosity models. They postulated shape factors for 1-, 2-, and 3D matrix blocks, as given in Table 1. In 1976, Kazemi et al. proposed different shape factors, which were derived using a finite-difference discretization. Kazemi et al. also postulated the generalization of the shape factor concept from single- to multiphase flow by introducing the phase relative permeability into Eq. 1. Thomas et al. (1983) found that they could accurately reproduce fine-grid single-porosity simulation results of water/oil countercurrent imbibition (in cubical blocks) if in their single-cell dual-porosity model they used a shape factor 25 / L2. In their dual-porosity simulation, however, they also used pseudorelative permeability curves and a pseudocapillary pressure, so it is not obvious whether the good fit was mainly caused by the shape factor they used, or by the pseudosaturation functions. Coats reported that the shape factor proposed by Kazemi is too low by a factor of 2, and derived new 1-, 2-, and 3D shape factors (1989); see Table 1. Ueda et al. (1989) also argued that the Kazemi shape factor should be multiplied by a factor 2 to 3, based on their work in which they compared dual porosity (two-phase) simulations with 1- and 2D fine-grid simulations. In 1993, Zimmerman et al. published a semi-analytical method for modeling of matrix-fracture flow in a dual-porosity model where the matrix blocks are modeled as spherical blocks (1993a). In their paper, they also show that the shape factor for spherical matrix blocks is given by p2 / R2 where R is the radius of the matrix block. In the same year, Chang derived an explicit formula for the single-phase shape factor for rectangular matrix blocks based on the full transient solution of the diffusion equation introducing new 1-, 2-, and 3D results to the shape-factor literature (1993). The same result was independently obtained in 1995 by Lim and Aziz. Both Chang and Lim and Aziz stressed that the shape factor, which had previously been regarded as a constant, is actually a function of time. In view of the wide spectrum of results and the apparent lack of consensus regarding which shape factor to use in simulations, a more detailed analysis into the reasons for the different shape factors cited in Table 1 seems desirable. We want to underline that in this paper we focus our attention to single-phase shape factors, thus avoiding the additional complications that arise in the discussion of two-phase matrix-fracture interaction because of relative permeability and capillary pressure. This allows us to more clearly illustrate the different approaches that the previously mentioned authors used.

Abstract Conventional hydrocarbon reservoirs, from an engineering and economic standpoint, are the easiest and most cost-efficient deposits to develop and produce. However, as economic deposits of conventional oil/gas become scarce, hydrocarbon recovered from tight sands and shale deposits will likely fill the void created by diminished conventional oil and gas sources. The purpose of this paper is to review the numerical methods available for simulating multiphase flow in highly fractured reservoirs and present a concise method to implement a fully implicit, two-phase numerical model for simulating multiphase flow, and predicting fluid recovery in highly fractured tight gas and shale gas reservoirs. The paper covers the five primary numerical modeling categories. It addresses the physical and theoretical concepts that support the development of numerical reservoir models and sequentially presents the stages of model development starting with mass balance fundamentals, Darcy’s law and the continuity equations. The paper shows how to develop and reduce the fluid transport equations. It also addresses equation discretization and linearization, model validation and typical model outputs. More advanced topics such as compositional models, reactive transport models, and artificial neural network models are also briefly discussed. The paper concludes with a discussion of field-scale model implementation challenges and constraints. The paper focuses on concisely and clearly presenting fundamental methods available to the novice petroleum engineer with the goal of improving their understanding of the inner workings of commercially available black box reservoir simulators. The paper assumes the reader has a working understanding of flow a porous media, Darcy’s law, and reservoir rock and fluid properties such as porosity, permeability, saturation, formation volume factor, viscosity, and capillary pressure. The paper does not explain these physical concepts neither are the laboratory tests needed to quantify these physical phenomena addressed. However, the paper briefly addresses these concepts in the context of sampling, uncertainty, upscaling, field-scale distribution, and the impact they have on field-scale numerical models.

Summary Foam injection is a promising enhanced-oil-recovery (EOR) technology that significantly improves the sweep efficiency of gas injection. Simulation of foam/oil displacement in reservoirs is an expensive process for conventional simulation because of the strongly nonlinear physics, such as multiphase flow and transport with oil/foam interactions. In this work, an operator-based linearization (OBL) approach, combined with the representation of foam by an implicit-texture (IT) model with two flow regimes, is extended for the simulation of the foam EOR process. The OBL approach improves the efficiency of the highly nonlinear foam-simulation problem by transforming the discretized nonlinear conservation equations into a quasilinear form using state-dependent operators. The state-dependent operators are approximated by discrete representation on a uniform mesh in parameter space. The numerical-simulation results are validated by using three-phasefractional-flow theory for foam/oil flow. Starting with an initial guess depending on the fitting of steady-state experimental data with oil, the OBL foam model is regressed to experimental observations using a gradient-optimization technique. A series of numerical validation studies is performed to investigate the accuracy of the proposed approach. The numerical model shows good agreement with analytical solutions at different conditions and with different foam parameters. With finer grids, the resolution of the simulation is better, but at the cost of more expensive computations. The foam-quality scan is accurately fitted to steady-state experimental data, except in the low-quality regime. In this regime, the used IT foam model cannot capture the upward-tilting pressure gradient (or apparent viscosity) contours. 1D and 3D simulation results clearly demonstrate two stages of foam propagation from inlet to outlet, as seen in the computed-tomography (CT) coreflood experiments: weak foam displaces most of the oil, followed by a propagation of stronger foam at lower oil saturation. OBL is a direct method to reduce nonlinearity in complex physical problems, which can significantly improve computational performance. Taking its accuracy and efficiency into account, the data-drivenOBL-based approach could serve as a platform for efficient numerical upscaling to field-scaleapplications.

AbstractThe Cobourg limestone is a heterogeneous argillaceous rock consisting of lighter nodular regions of calcite and dolomite, interspersed with darker regions composed of calcite, dolomite, quartz and a clay fraction. The intact permeability of the Cobourg limestone is estimated to be in the range of K ∈ (10−23, 10−19) m2. This paper discusses the factors influencing the measurement of the intact permeability of the Cobourg limestone and presents an upscaling approach for estimating this parameter. The procedure first involves the dissection of a cuboidal sample of the rock measuring, 80 mm × 120 mm × 300 mm, into ten 8 mm-thick slabs. Digital imaging and mapping of the larger surfaces of these sections are used to create, from both surface image extrusion and surface image interpolation techniques, the fabric within the dissected regions. The estimated permeabilities of the lighter and darker regions are used in the computational models of the computer-generated fabric to estimate the effective permeability of the rock. These results are complemented by estimates derived from mathematical theories for estimating permeabilities of multiphasic composites.