Multiscale mixed/mimetic methods on corner-point grids

Multiscale simulation is a promising approach to facilitate direct simulation of large and complex grid models for highly heterogeneous petroleum reservoirs. Unlike traditional simulation, approaches based on upscaling/downscaling, multiscale methods seek to solve the full flow problem by incorporating subscale heterogeneities into local discrete approximation spaces. We consider a multiscale formulation based on a hierarchical grid approach, where basis functions with subgrid resolution are computed numerically to correctly and accurately account for subscale variations from an underlying (fine-scale) geomodel when solving the global flow equations on a coarse grid. By using multiscale basis functions to discretise the global flow equations on a (moderately sized) coarse grid, one can retain the efficiency of an upscaling method and, at the same time, produce detailed and conservative velocity fields on the underlying fine grid. For pressure equations, the multiscale mixed finite-element method (MsMFEM) has been shown to be a particularly versatile approach. In this paper, we extend the method to corner-point grids, which is the industry standard for modelling complex reservoir geology. To implement MsMFEM, one needs a discretisation method for solving local flow problems on the underlying fine grids. In principle, any stable and conservative method can be used. Here, we use a mimetic discretisation, which is a generalisation of mixed finite elements that gives a discrete inner product, allows for polyhedral elements, and can (easily) be extended to curved grid faces. The coarse grid can, in principle, be any partition of the subgrid, where each coarse block is a connected collection of subgrid cells. However, we argue that, when generating coarse grids, one should follow certain simple guidelines to achieve improved accuracy. We discuss partitioning in both index space and physical space and suggest simple processing techniques. The versatility and accuracy of the new multiscale mixed methodology is demonstrated on two corner-point models: a small Y-shaped sector model and a complex model of a layered sedimentary bed. A variety of coarse grids, both violating and obeying the above mentioned guidelines, are employed. The MsMFEM solutions are compared with a reference solution obtained by direct simulation on the subgrid.     

Integration of local–global upscaling and grid adaptivity for simulation of subsurface flow in heterogeneous formations

We propose a methodology, called multilevel local–global (MLLG) upscaling, for generating accurate upscaled models of permeabilities or transmissibilities for flow simulation on adapted grids in heterogeneous subsurface formations. The method generates an initial adapted grid based on the given fine-scale reservoir heterogeneity and potential flow paths. It then applies local–global (LG) upscaling for permeability or transmissibility [7], along with adaptivity, in an iterative manner. In each iteration of MLLG, the grid can be adapted where needed to reduce flow solver and upscaling errors. The adaptivity is controlled with a flow-based indicator. The iterative process is continued until consistency between the global solve on the adapted grid and the local solves is obtained. While each application of LG upscaling is also an iterative process, this inner iteration generally takes only one or two iterations to converge. Furthermore, the number of outer iterations is bounded above, and hence, the computational costs of this approach are low. We design a new flow-based weighting of transmissibility values in LG upscaling that significantly improves the accuracy of LG and MLLG over traditional local transmissibility calculations. For highly heterogeneous (e.g., channelized) systems, the integration of grid adaptivity and LG upscaling is shown to consistently provide more accurate coarse-scale models for global flow, relative to reference fine-scale results, than do existing upscaling techniques applied to uniform grids of similar densities. Another attractive property of the integration of upscaling and adaptivity is that process dependency is strongly reduced, that is, the approach computes accurate global flow results also for flows driven by boundary conditions different from the generic boundary conditions used to compute the upscaled parameters. The method is demonstrated on Cartesian cell-based anisotropic refinement (CCAR) grids, but it can be applied to other adaptation strategies for structured grids and extended to unstructured grids.     

Grid adaptation for the Dirichlet–Neumann representation method and the multiscale mixed finite-element method

A Dirichlet–Neumann representation method was recently proposed for upscaling and simulating flow in reservoirs. The DNR method expresses coarse fluxes as linear functions of multiple pressure values along the boundary and at the center of each coarse block. The number of flux and pressure values at the boundary can be adjusted to improve the accuracy of simulation results and, in particular, to resolve important fine-scale details. Improvement over existing approaches is substantial especially for reservoirs that contain high-permeability streaks or channels. As an alternative, the multiscale mixed finite-element (MsMFE) method was designed to obtain fine-scale fluxes at the cost of solving a coarsened problem, but can also be used as upscaling methods that are flexible with respect to geometry and topology of the coarsened grid. Both methods can be expressed in mixed-hybrid form, with local stiffness matrices obtained as “inner products” of numerically computed basis functions with fine-scale sub-resolution. These basis functions are determined by solving local flow problems with piecewise linear Dirichlet boundary conditions for the DNR method and piecewise constant Neumann conditions for MsMFE. Adding discrete pressure points in the DNR method corresponds to subdividing faces in the coarse grid and hence increasing the number of basis functions in the MsMFE method. The methods show similar accuracy for 2D Cartesian cases, but the MsMFE method is more straightforward to formulate in 3D and implement for general grids.     

A comparison of multiscale methods for elliptic problems in porous media flow

We review and perform comparison studies for three recent multiscale methods for solving elliptic problems in porous media flow; the multiscale mixed finite-element method, the numerical subgrid upscaling method, and the multiscale finite-volume method. These methods are based on a hierarchical strategy, where the global flow equations are solved on a coarsened mesh only. However, for each method, the discrete formulation of the partial differential equations on the coarse mesh is designed in a particular fashion to account for the impact of heterogeneous subgrid structures of the porous medium. The three multiscale methods produce solutions that are mass conservative on the underlying fine mesh. The methods may therefore be viewed as efficient, approximate fine-scale solvers, i.e., as an inexpensive alternative to solving the elliptic problem on the fine mesh. In addition, the methods may be utilized as an alternative to upscaling, as they generate mass-conservative solutions on the coarse mesh. We therefore choose to also compare the multiscale methods with a state-of-the-art upscaling method – the adaptive local–global upscaling method, which may be viewed as a multiscale method when coupled with a mass-conservative downscaling procedure. We investigate the properties of all four methods through a series of numerical experiments designed to reveal differences with regard to accuracy and robustness. The numerical experiments reveal particular problems with some of the methods, and these will be discussed in detail along with possible solutions. Next, we comment on implementational aspects and perform a simple analysis and comparison of the computational costs associated with each of the methods. Finally, we apply the three multiscale methods to a dynamic two-phase flow case and demonstrate that high efficiency and accurate results can be obtained when the subgrid computations are made part of a preprocessing step and not updated, or updated infrequently, throughout the simulation. The research is funded by the Research Council of Norway under grant nos. 152732 and 158908.     

Upscaling Hydraulic Conductivities in Cross-Bedded Formations

Modern geostatistical techniques allow the generation of high-resolution heterogeneous models of hydraulic conductivity containing millions to billions of cells. Selective upscaling is a numerical approach for the change of scale of fine-scale hydraulic conductivity models into coarser scale models that are suitable for numerical simulations of groundwater flow and mass transport. Selective upscaling uses an elastic gridding technique to selectively determine the geometry of the coarse grid by an iterative procedure. The geometry of the coarse grid is built so that the variances of flow velocities within the coarse blocks are minimum. Selective upscaling is able to handle complex geological formations and flow patterns, and provides full hydraulic conductivity tensor for each block. Selective upscaling is applied to a cross-bedded formation in which the fine-scale hydraulic conductivities are full tensors with principal directions not parallel to the statistical anisotropy of their spatial distribution. Mass transport results from three coarse-scale models constructed by different upscaling techniques are compared to the fine-scale results for different flow conditions. Selective upscaling provides coarse grids in which mass transport simulation is in good agreement with the fine-scale simulations, and consistently superior to simulations on traditional regular (equal-sized) grids or elastic grids built without accounting for flow velocities.     

Multiscale mass conservative domain decomposition preconditioners for elliptic problems on irregular grids

Multiscale methods can in many cases be viewed as special types of domain decomposition preconditioners. The localisation approximations introduced within the multiscale framework are dependent upon both the heterogeneity of the reservoir and the structure of the computational grid. While previous works on multiscale control volume methods have focused on heterogeneous elliptic problems on regular Cartesian grids, we have tested the multiscale control volume formulations on two-dimensional elliptic problems involving heterogeneous media and irregular grid structures. Our study shows that the tangential flow approximation commonly used within multiscale methods is not suited for problems involving rough grids. We present a more robust mass conservative domain decomposition preconditioner for simulating flow in heterogeneous porous media on general grids.     

Upscaling facies models to preserve connectivity of designated facies

An upscaling algorithm has been developed that generates an irregular coarse grid that preserves flow connectivity by applying a rule-based upscaling algorithm to a fine-scale facies distribution. The algorithm is demonstrated using stochastically generated paleo-fluvial facies distributions. First, an irregular grid honoring the channel facies is created, followed by computation of effective anisotropic parameters for all coarse-grid cells. For the apparent layer-cake geometry of overbank deposits seen in outcrop, two local upscaling methods are compared: (1) the layered system approximation and (2) the mode. To assess upscaling performance, flow simulations for the original and upscaled grids are compared. The horizontal layered approximation (arithmetic mean) performs poorly, over-predicting lateral connectivity where even infrequent disconnection becomes important. Performance of the mode as an upscaling algorithm depends on the probability that a coarse-grid cell will be dominated by a single facies, and it performs surprisingly well because the upscaled grid-generation algorithm honors the channels, informing the upscaling process. Lastly, the irregular coarse grid was compared to a uniform coarse grid, showing superior performance with the irregular grid. The reduction in grid size achieved by irregular-grid generation will be a function of the geometrical complexity of the geologic objects to be honored.     

A multiscale method for distributed parameter estimation with application to reservoir history matching

A method for multiscale parameter estimation with application to reservoir history matching is presented. Starting from a given fine-scale model, coarser models are generated using a global upscaling technique where the coarse models are tuned to match the solution of the fine model. Conditioning to dynamic data is done by history-matching the coarse model. Using consistently the same resolution both for the forward and inverse problems, this model is successively refined using a combination of downscaling and history matching until model-matching dynamic data are obtained at the finest scale. Large-scale corrections are obtained using fast models, which, combined with a downscaling procedure, provide a better initial model for the final adjustment on the fine scale. The result is thus a series of models with different resolution, all matching history as good as possible with this grid. Numerical examples show that this method may significantly reduce the computational effort and/or improve the quality of the solution when achieving a fine-scale match as compared to history-matching directly on the fine scale.     

Statistical assignment of upscaled flow functions for an ensemble of geological models

Upscaled flow functions are often needed to account for the effects of fine-scale permeability heterogeneity in coarse-scale simulation models. We present procedures in which the required coarse-scale flow functions are statistically assigned to an ensemble of upscaled geological models. This can be viewed as an extension and further development of a recently developed ensemble level upscaling (EnLU) approach. The method aims to efficiently generate coarse-scale flow models capable of reproducing the ensemble statistics (e.g., cumulative distribution function) of fine-scale flow predictions for multiple reservoir models. The most expensive part of standard coarsening procedures is typically the generation of upscaled two-phase flow functions (e.g., relative permeabilities). EnLU provides a means for efficiently generating these upscaled functions using stochastic simulation. This involves the use of coarse-block attributes that are both fast to compute and correlate closely with the upscaled two-phase functions. In this paper, improved attributes for use in EnLU, namely the coefficient of variation of the fine-scale single-phase velocity field (computed during computation of upscaled absolute permeability) and the integral range of the fine-scale permeability variogram, are identified. Geostatistical simulation methods, which account for spatial correlations of the statistically generated upscaled functions, are also applied. The overall methodology thus enables the efficient generation of coarse-scale flow models. The procedure is tested on 3D well-driven flow problems with different permeability distributions and variable fluid mobility ratios. EnLU is shown to capture the ensemble statistics of fine-scale flow results (water and oil flow rates as a function of time) with similar accuracy to full flow-based upscaling methods but with computational speedups of more than an order of magnitude.     

Nodal and Mixed Finite Elements for the Numerical Homogenization of 3D Permeability

The aim of upscaling is to determine equivalent homogeneous parameters at a coarse-scale from a spatially oscillating fine-scale parameter distribution. To be able to use a limited number of relatively large grid-blocks in numerical oil reservoir simulators or groundwater models, upscaling of the permeability is frequently applied. The spatial fine-scale permeability distribution is generally obtained from geological and geostatistical models. After upscaling, the coarse-scale permeabilities are incorporated in the relatively large grid-blocks of the numerical model. If the porous rock may be approximated as a periodic medium, upscaling can be performed by the method of homogenization. In this paper the homogenization is performed numerically, which gives rise to an approximation error. The complementarity between two different numerical methods – the conformal-nodal finite element method and the mixed-hybrid finite element method – has been used to quantify this error. These two methods yield respectively upper and lower bounds for the eigenvalues of the coarse-scale permeability tensor. Results of 3D numerical experiments are shown, both for the far field and around wells.     

Mixed multiscale finite element methods using approximate global information based on partial upscaling

The use of limited global information in multiscale simulations is needed when there is no scale separation. Previous approaches entail fine-scale simulations in the computation of the global information. The computation of the global information is expensive. In this paper, we propose the use of approximate global information based on partial upscaling. A requirement for partial homogenization is to capture long-range (non-local) effects present in the fine-scale solution, while homogenizing some of the smallest scales. The local information at these smallest scales is captured in the computation of basis functions. Thus, the proposed approach allows us to avoid the computations at the scales that can be homogenized. This results in coarser problems for the computation of global fields. We analyze the convergence of the proposed method. Mathematical formalism is introduced, which allows estimating the errors due to small scales that are homogenized. The proposed method is applied to simulate two-phase flows in heterogeneous porous media. Numerical results are presented for various permeability fields, including those generated using two-point correlation functions and channelized permeability fields from the SPE Comparative Project (Christie and Blunt, SPE Reserv Evalu Eng 4:308–317, 2001). We consider simple cases where one can identify the scales that can be homogenized. For more general cases, we suggest the use of upscaling on the coarse grid with the size smaller than the target coarse grid where multiscale basis functions are constructed. This intermediate coarse grid renders a partially upscaled solution that contains essential non-local information. Numerical examples demonstrate that the use of approximate global information provides better accuracy than purely local multiscale methods.     

Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible,three-phase flow with gravity

Most practical reservoir simulation studies are performed using the so-called black oil model, in which the phase behavior is represented using solubilities and formation volume factors. We extend the multiscale finite-volume (MSFV) method to deal with nonlinear immiscible three-phase compressible flow in the presence of gravity and capillary forces (i.e., black oil model). Consistent with the MSFV framework, flow and transport are treated separately and differently using a sequential implicit algorithm. A multiscale operator splitting strategy is used to solve the overall mass balance (i.e., the pressure equation). The black-oil pressure equation, which is nonlinear and parabolic, is decomposed into three parts. The first is a homo geneous elliptic equation, for which the original MSFV method is used to compute the dual basis functions and the coarse-scale transmissibilities. The second equation accounts for gravity and capillary effects; the third equation accounts for mass accumulation and sources/ sinks (wells). With the basis functions of the elliptic part, the coarse-scale operator can be assembled. The gravity/capillary pressure part is made up of an elliptic part and a correction term, which is computed using solutions of gravity-driven local problems. A particular solution represents accumulation and wells. The reconstructed fine-scale pressure is used to compute the fine-scale phase fluxes, which are then used to solve the nonlinear saturation equations. For this purpose, a Schwarz iterative scheme is used on the primal coarse grid. The framework is demonstrated using challenging black-oil examples of nonlinear compressible multiphase flow in strongly heterogeneous formations.     

A system-theoretical approach to selective grid coarsening of reservoir models

From a system-theoretical point of view and for a given configuration of wells, there are only a limited number of degrees of freedom in the input–output dynamics of a reservoir system. This means that a large number of combinations of the state variables (pressure and saturation values) are not actually controllable and observable from the wells, and accordingly, they are not affecting the input–output behavior of the system. In an earlier publication, we therefore proposed a control-relevant upscaling methodology that uniformly coarsens the reservoir. Here, we present a control-relevant selective (i.e. non-uniform) coarsening (CRSC) method, in which the criterion for grid size adaptation is based on ranking the grid block contributions to the controllability and observability of the reservoir system. This multi-level CRSC method is attractive for use in iterative procedures such as computer-assisted flooding optimization for a given configuration of wells. In contrast to conventional flow-based coarsening techniques our method is independent of the specific flow rates or pressures imposed at the wells. Moreover the system-theoretical norms employed in our method provide tight upper bounds to the ‘input–output energy’ of the fine and coarse systems. These can be used as an a priori error-estimate of the performance of the coarse model. We applied our algorithm to two numerical examples and found that it can accurately reproduce results from the corresponding fine-scale simulations, while significantly speeding up the simulation.     

FV-MHMM method for reservoir modeling

The present paper proposes a new family of multiscale finite volume methods. These methods usually deal with a dual mesh resolution, where the pressure field is solved on a coarse mesh, while the saturation fields, which may have discontinuities, are solved on a finer reservoir grid, on which petrophysical heterogeneities are defined. Unfortunately, the efficiency of dual mesh methods is strongly related to the definition of up-gridding and down-gridding steps, allowing defining accurately pressure and saturation fields on both fine and coarse meshes and the ability of the approach to be parallelized. In the new dual mesh formulation we developed, the pressure is solved on a coarse grid using a new hybrid formulation of the parabolic problem. This type of multiscale method for pressure equation called multiscale hybrid-mixed method (MHMM) has been recently proposed for finite elements and mixed-finite element approach (Harder et al. 2013). We extend here the MH-mixed method to a finite volume discretization, in order to deal with large multiphase reservoir models. The pressure solution is obtained by solving a hybrid form of the pressure problem on the coarse mesh, for which unknowns are fluxes defined on the coarse mesh faces. Basis flux functions are defined through the resolution of a local finite volume problem, which accounts for local heterogeneity, whereas pressure continuity between cells is weakly imposed through flux basis functions, regarded as Lagrange multipliers. Such an approach is conservative both on the coarse and local scales and can be easily parallelized, which is an advantage compared to other existing finite volume multiscale approaches. It has also a high flexibility to refine the coarse discretization just by refinement of the lagrange multiplier space defined on the coarse faces without changing nor the coarse nor the fine meshes. This refinement can also be done adaptively w.r.t. a posteriori error estimators. The method is applied to single phase (well-testing) and multiphase flow in heterogeneous porous media.     

A scalable multistage linear solver for reservoir models with multisegment wells

Wellbore flow and interactions between wells and the reservoir can be complex. Accurate modeling of these behaviors is especially important for multilateral and other advanced wells. This paper describes a new scalable linear solver for flow simulation of detailed reservoir models with advanced wells and well groups. A general purpose research simulator serves as the computational platform, in which a multisegment well (MsWell) model is used to describe wellbore flow. In the MsWell model, the wellbore is discretized into a number of segments. Hence, the MsWell model adds a large number of equations and unknowns, which are fully coupled to the reservoir model. Operating constraints on groups of wells add one more level of complexity to the system. The new linear solver is a generalized two-stage constrained pressure residual preconditioner. A global pressure system is obtained algebraically in the first stage. The system represents the pressure coupling between the reservoir and wells accurately. The well groups are disaggregated into individual multisegment wells, which then are further reduced to a standard well-like form. The two-stage scheme serves as the inner loop of a generalized minimum residual solver. Algebraic multigrid is used to compute the first-stage pressure solution; a special block-based incomplete lower–upper preconditioner is used for the second stage. We demonstrate the superior performance of this new solver compared with state-of-the-art methods using a variety of highly detailed reservoir models with complex wells and well groups.     

Finite element modelling of three-phase flow in deforming saturated oil reservoirs

A numerical simulation is presented for three-dimensional three-phase fluid flow in a deforming saturated oil reservoir. The mathematical formulation describes a fully coupled governing equation systen which consists of the equilibrium and continuity equations for three immiscible fluids flowing in a porous medium. An elastoplastic soil model, based on a Mohr–Coulomb yield surface, is used. The finite element method is applied to obtain simultaneous solutions to the governing equations where displacement and fluid pressures are the primary unknowns. The final discretized equations are solved by a direct solver using fully implicit procedures. The developed model is used to investigate the problems of three-phase fluid flow in a deforming saturated oil reservoir.     

An upscaling approach using adaptive multi-resolution upgridding and automated relative permeability adjustment

The upscaling process of a high-resolution geostatistical reservoir model to a dynamic simulation grid model plays an important role in a reservoir study. Several upscaling methods have been proposed in order to create balance between the result accuracy and computation speed. Usually, a high-resolution grid model is upscaled according to the heterogeneities assuming single phase flow. However, during injection processes, the relative permeability adjustment is required. The so-called pseudo-relative permeability curves are accepted, if their corresponding coarse model is a good representation of the fine-grid model. In this study, an upscaling method based on discrete wavelet transform (WT) is developed for single-phase upscaling based on the multi-resolution analysis (MRA) concepts. Afterwards, an automated optimization method is used in which evolutionary genetic algorithm is applied to estimate the pseudo-relative permeability curves described with B-spline formulation. In this regard, the formulation of B-spline is modified in order to describe the relative permeability curves. The proposed procedure is evaluated in the gas injection case study from the SPE 10th comparative solution project’s data set which provides a benchmark for upscaling problems [1]. The comparisons of the wavelet-based upscaled model to the high-resolution model and uniformly coarsened model show considerable speedup relative to the fine-grid model and better accuracy relative to the uniformly coarsened model. In addition, the run time of the wavelet-based coarsened model is comparable with the run time of the uniformly upscaled model. The optimized coarse models increase the speed of simulation up to 90% while presenting similar results as fine-grid models. Besides, using two different production/injection scenarios, the superiority of WT upscaling plus relative permeability adjustment over uniform upscaling and relative permeability adjustment is presented. This study demonstrates the proposed upscaling workflow as an effective tool for a reservoir simulation study to reduce the required computational time.     

Near-well upscaling for three-phase flows

Large-scale flow models constructed using standard coarsening procedures may not accurately resolve detailed near-well effects. Such effects are often important to capture, however, as the interaction of the well with the formation can have a dominant impact on process performance. In this work, a near-well upscaling procedure, which provides three-phase well-block properties, is developed and tested. The overall approach represents an extension of a recently developed oil–gas upscaling procedure and entails the use of local well computations (over a region referred to as the local well model (LWM)) along with a gradient-based optimization procedure to minimize the mismatch between fine and coarse-scale well rates, for oil, gas, and water, over the LWM. The gradients required for the minimization are computed efficiently through solution of adjoint equations. The LWM boundary conditions are determined using an iterative local-global procedure. With this approach, pressures and saturations computed during a global coarse-scale simulation are interpolated onto LWM boundaries and then used as boundary conditions for the fine-scale LWM computations. In addition to extending the overall approach to the three-phase case, this work also introduces new treatments that provide improved accuracy in cases with significant flux from the gas cap into the well block. The near-well multiphase upscaling method is applied to heterogeneous reservoir models, with production from vertical and horizontal wells. Simulation results illustrate that the method is able to accurately capture key near-well effects and to provide predictions for component production rates that are in close agreement with reference fine-scale results. The level of accuracy of the procedure is shown to be significantly higher than that of a standard approach which uses only upscaled single-phase flow parameters.     

A multi-resolution workflow to generate high-resolution models constrained to dynamic data

Distance-based stochastic techniques have recently emerged in the context of ensemble modeling, in particular for history matching, model selection and uncertainty quantification. Starting with an initial ensemble of realizations, a distance between any two models is defined. This distance is defined such that the objective of the study is incorporated into the geological modeling process, thereby potentially enhancing the efficacy of the overall workflow. If the intent is to create new models that are constrained to dynamic data (history matching), the calculation of the distance requires flow simulation for each model in the initial ensemble. This can be very time consuming, especially for high-resolution models. In this paper, we present a multi-resolution framework for ensemble modeling. A distance-based procedure is employed, with emphasis on the rapid construction of multiple models that have improved dynamic data conditioning. Our intent is to construct new high-resolution models constrained to dynamic data, while performing most of the flow simulations only on upscaled models. An error modeling procedure is introduced into the distance calculations to account for potential errors in the upscaling. Based on a few fine-scale flow simulations, the upscaling error is estimated for each model using a clustering technique. We demonstrate the efficiency of the method on two examples, one where the upscaling error is small, and another where the upscaling error is significant. Results show that the error modeling procedure can accurately capture the error in upscaling, and can thus reproduce the fine-scale flow behavior from coarse-scale simulations with sufficient accuracy (in terms of uncertainty predictions). As a consequence, an ensemble of high-resolution models, which are constrained to dynamic data, can be obtained, but with a minimum of flow simulations at the fine scale.     

Numerical Homogenization of Two-Phase Flow in Porous Media

Homogenization has proved its effectiveness as a method of upscaling for linear problems, as they occur in single-phase porous media flow for arbitrary heterogeneous rocks. Here we extend the classical homogenization approach to nonlinear problems by considering incompressible, immiscible two-phase porous media flow. The extensions have been based on the principle of preservation of form, stating that the mathematical form of the fine-scale equations should be preserved as much as possible on the coarse scale. This principle leads to the required extensions, while making the physics underlying homogenization transparent. The method is process-independent in a way that coarse-scale results obtained for a particular reservoir can be used in any simulation, irrespective of the scenario that is simulated. Homogenization is based on steady-state flow equations with periodic boundary conditions for the capillary pressure. The resulting equations are solved numerically by two complementary finite element methods. This makes it possible to assess a posteriori error bounds.