No-flow Boundary Conditions Research Articles

1-D steady state runoff production in light of queuing theory: Heterogeneity, connectivity, and scale

[1] We used the frameworks of queuing theory and connectivity to study the runoff generated under constant rainfall on a one-dimensional slope with randomly distributed infiltrability. The equivalence between the stationary runoff-runon equation and the customers waiting time in a single server queue provides a theoretical link between the statistical description of infiltrability and that of runoff flow rate. Five distributions of infiltrability, representing soil heterogeneities at different scales, are considered: four uncorrelated (exponential, bimodal, lognormal, uniform) and one autocorrelated (lognormal, with or without a nugget). The existing theoretical results are adapted to the hydrological framework for the exponential case, and new theoretical developments are proposed for the bimodal law. Numerical simulations validate these results and improve our understanding of runoff-runon for all of the distributions. The quantities describing runoff generation (runoff one-point statistics) and its organization into patterns (patterns statistics and connectivity) are studied as functions of rainfall rate. The variables describing the wet areas are also compared to those describing the rainfall excess areas, i.e., the areas where rainfall exceeds infiltrability. Preliminary results concerning the structural and functional connectivity functions are provided, as well as a discussion about the origin of scale effects in such a system. We suggest that the upslope no-flow boundary condition may be responsible for the dependence of the runoff coefficient on the scale of observation. Queuing theory appears to be a promising framework for runoff-runon modeling and hydrological connectivity problems.

Superposition of steady-state coning solutions for multiple wells with reservoir boundaries

The critical oil (or gas) rate to avoid coning of unwanted fluids into production wells is an important design parameter. Simulation methods are typically used to predict critical rates in reservoirs with complex heterogeneities and boundaries, but they are manpower intensive and prone to errors when large grid blocks are used. Current analytical methods are useful to provide benchmark solutions for simulation, but their assumptions are too restrictive for more general use. Thus, there is a need for improved analytical methods to account for well patterns and more complex boundaries that can serve as additional benchmarks for simulation, and in some cases to predict critical rates. This paper makes analytical solutions more realistic by extending existing single-well analytical solutions to account for multiple wells and common no-flow and constant pressure boundary conditions. A potential function is derived to superimpose existing single well coning solutions for simultaneous two-phase flow. Capillary pressure and relative permeability effects on coning are included in the new potential function. Vertical equilibrium (VE) and steady-state flow are assumed. Comparisons with finely-gridded simulation show good agreement in predicted critical oil rates when steady state and VE are approached. VE and steady-state are approximately achieved when aspect ratios are greater than about 10. The predicted critical rates are conservative for aspect ratios less than 10. The new solutions demonstrate for several examples how the critical rates from one well are affected by reservoir boundaries, and the two-phase production rates of other wells. • Developed new analytical solutions for coning when two phases are flowing simultaneously. • Developed a potential function that allows for superposition so that multiple wells and no-flow and constant potential reservoir boundary conditions can be analyzed. • Demonstrates that boundary conditions and other wells can significantly alter the critical oil rates to avoid coning.

Conditions for lateral downslope unsaturated flow and effects of slope angle on soil moisture movement

► Lateral downslope flow at slope surface and bottom is due to the state of drying. ► Speed of lateral flow to become parallel to slope surface is due to drainage rate. ► Lateral flow component parallel to slope surface is more sensitive to slope angle. ► Influence of the slope angle is the most significant in the middle soil layer. The necessary conditions for lateral downslope flow in the unsaturated zone in hillslopes have not been clearly identified, and the effects of different slope angles on soil moisture movement are rarely studied. These questions were investigated through drying processes in a homogeneous and isotropic sloping soil tank in this study. Three experiments with different slope angles were conducted. The results showed that the flow direction in the unsaturated zone gradually rotates counterclockwise from the vertical direction to the lateral downslope direction parallel to the slope surface during the drying phase. The lateral downslope flow parallel to the slope surface first appears at the surface (0–5 cm) and the bottom of the sandy loam in the tank. The lateral downslope flow parallel to the slope surface at the surface is due to the state of drying rather than the no-flow upper boundary condition. The lateral downslope flow parallel to the slope surface at the bottom of the sandy loam is also caused by the state of drying. More importantly, the results of this study showed that the rate at which the lateral downslope flow turns to be parallel to the slope surface is affected by the drainage rate. In contrast to the lateral downslope flow component parallel to the slope surface, the flow component normal to the slope surface is less sensitive to changes in the slope angle. The influence of the slope angle on the flow component normal to the slope surface is greatest in the middle layer.

Upscaling of upward CO2 migration in 2D system

A procedure for upscaling CO2 buoyancy driven upward migration in finite-difference simulation models is presented in this work. This upscaling procedure accounts for capillary and buoyancy forces to enable CO2 upward migration modeling in coarser grids while accounting for dominant fine-scaled geological effects. The developed method is applied to 2D domains with no-flow boundary conditions. The absolute permeability field is correlated in the horizontal direction, with zero correlation in the vertical direction. Capillary pressure is parameterized using a Leveret J-function. A Dykstra-Parsons coefficient of 0.7 was used to generate a relatively heterogeneous absolute permeability field and hence test the developed algorithm under more stringent conditions. Multiphase flow upscaling is improved by accounting for spatial connectivity (percolation), which enables us to obtain more realistic rock-fluid pseudo-functions and capture effects of local capillary trapping at the fine scale (meso-scale trapping). The upscaling method and estimation of rock-fluid functions are numerically tested and compared with currently accepted single and multiphase flow upscaling methods. Results show that single-phase flow upscaling is insufficient, because it fails to adequately predict mobility and residual saturation, and hence multiphase flow upscaling should be employed. Significant improvement in gas travel time (representative of mobility) and trapped CO2 saturation (representative of trapped saturation) are observed when spatial connectivity (percolation) is included. The simulation execution time reduces 17-fold through upscaling. This speedup will enable simulating 3D CO2 sequestration simulation scenarios.

Note on the uniqueness of subsonic Euler flows in an axisymmetric nozzle

In this paper, we consider the uniqueness of globally subsonic compressible flows through an infinitely long axisymmetric nozzle. The flow is governed by the steady Euler equations and satisfies no-flow boundary conditions on the nozzle walls. We will show that for given mass flux and Bernoulli’s function in the upstream, the subsonic flow is unique in the class of all axisymmetric solutions, which possess the asymptotic behaviors at the far fields. This result extends the uniqueness of solutions in the previous paper Du and Duan (2011) [1].

Time-Lapse Variations of Multi-Component Electrical Resistivity Measurements Acquired In High-Angle Wells

Measurements of electrical conductivity anisotropy are becoming increasingly common in LWD and open-hole environments. At the same time, field experience indicates that invasion and time of logging can have a significant effect on multi-component induction measurements acquired in high-angle wells, specifically in the presence of highly permeable formations and large fluid density and viscosity contrasts. Some efforts have been made to study these effects, but the actual spatial distribution of fluid saturation due to invasion has not been considered previously. We use a commercial multiphase fluid-flow simulator to reproduce invasion behavior in high-angle wells. Simulations provide an effective and reliable environment to study time-lapse effects of different geometrical and petrophysical parameters on the ensuing spatial distribution of fluid saturation. The study considers the effects of both water- and oil-base muds invading hydrocarbon-bearing formations in interbedded sand-shale sequences. Simulations show that gravity segregation causes off-centered and asymmetric spatial distributions of electrical resistivity in high-angle wells. Further, no-flow boundary conditions at impermeable shale boundaries cause high local concentrations of mud-filtrate. These effects can vary significantly with time. Depending on the time of measurement acquisition, three-dimensional (3D) effects due to invasion can distort multi-component apparent resistivity estimates. Consequently, inversion techniques should not be blindly used to synthesize all the measurements into estimates of conductivity anisotropy. It is necessary that 3D effects, including those due to invasion, be diagnosed in the analysis prior to estimating anisotropic conductivity parallel and perpendicular to bedding plane. Our work suggests that cross-coupled components of opposite sign and non-monotonic time variations in the coplanar components could diagnose time dependent invasion effects on multi-component induction measurements.

Analytic solutions with ionic flow for a pressure transmission test on shale

Ion/water transport between drilling fluid and shale formations can result in altered pore pressure in the shale. This chemically induced pore pressure plays an important role in controlling shale stability during drilling operations. In this paper analytic solutions are derived for a transient pressure transmission test on shale by solving the previously derived partially-coupled diffusivity equations taking into account chemical ionic flow effects. The diffusivity equations were solved under a no-flow boundary condition which is prescribed in a particular kind of transient pressure transmission test. By determining the three chemo-poroelastic coefficients K I (permeability coefficient), K II (membrane efficiency coefficient), and D eff (ion diffusivity) experimentally, the analytic solutions are compared with the actual experimental data. In general, a good agreement between the analytic solutions and the actual pressure measurements in the time-dependent pressure transmission tests has been found and is presented. The analytic solutions are provided for both the linear and cylindrical pressure transmission geometries. The pressure transmission experimental data discussed in this paper are from linear-geometry tests. Results presented in this paper can be used to study the transient pore pressure and stresses in the near-wellbore region, and to optimize the drilling fluid design (e.g., optimized salt type, salt concentration, optimized membrane efficiency, fluid pressure etc.).

Universal cokriging of hydraulic heads accounting for boundary conditions

When contouring scalar potentials from point observations the process can often benefit from including the known effects of boundary curves with specified potential or gradient. Here we consider the hydraulic head in an aquifer and both no-flow and constant-head boundary conditions. We present a new approach to enforcing that equipotential contours be normal to no-flow boundaries. A constant-head boundary, with unknown head, can be included through the same process by rotating the boundary vector by 90°. Collocated observations of heads and boundaries can specify a constant-head boundary of known value. We estimate head given both head and boundary condition observations, cokriging with both types of information. Our new approach uses gradient vectors in contrast with previous approximate finite-difference methods that include boundary conditions in kriging . Either the approach given here or the finite-difference method must be implemented with smooth covariance models, e.g., Gaussian, generalized Cauchy, and Matérn.

Numerical simulation of the effect of the sloping submarine outlet-capping on tidal groundwater head fluctuation in confined coastal aquifers

The submarine outlet of a coastal confined aquifer is usually covered by less-permeable material such as silt and fine sand which forms a sloping or almost flat seabed. Previous analytical studies of the tide-induced head fluctuation usually assumed a vertical interface between the outlet-capping and the seawater and neglected the effects of the outlet-capping's slope. Here we conducted a series of numerical simulations to investigate the effect of the slope of the outlet-capping on the tide-induced head fluctuations in a coastal confined aquifer. The numerical simulations demonstrated that when the hydraulic diffusivity, slope and/or the outlet-capping leakance are large, the tidal loading effect is relatively weak and the leakance dominates. In this case the analytical solution applies. In general, the joint actions of the outlet-capping leakance and the tidal loading effects result in complicated, 2-dimensional flow in the aquifer near the shoreline. For aquifers with small hydraulic diffusivity, small slope and/or small outlet-capping leakance, the outlet acts approximately as a no-flow boundary condition and the tidal loading dominates. In this case, negative phase shifts (time-advance) may occur near the coastline. The study provides potential guide to infer the aquifer's submarine structure and parameter from tidal head fluctuation observations in inland wells.

One-, two-, and three-dimensional solutions for counter-current steady-state two-phase subsurface flow

This paper derives analytical solutions for steady-state one-dimensional (1-D), two-dimensional (2-D), and three-dimensional (3-D) two-phase immiscible subsurface flow for a counter-current problem. Since the governing equations are highly nonlinear, 2-D and 3-D derivations are generally difficult to obtain. The primary purpose for the solutions is to test finite difference/volume/element computer programs for accuracy and scalability using architectures ranging from PCs to parallel high performance computers. This derivation is accomplished by first solving for the saturation of water in terms of a function that is a solution to Laplace’s equation to achieve a set of partial differential equations that allows some degree of latitude in the choice of boundary conditions. Separation of variables and Fourier series are used to obtain the final solution. The test problem consists of a rectangular block of soil where specified pressure is applied at the top and bottom of the sample, and no-flow boundary conditions are imposed on the sides. The pressure at the top of the sample is a step function that allows the testing of adaptive meshing or concentration of grid points in action zones.

Barotropic Instability of a Zonal Jet: From Nondivergent Perturbations on the β Plane to Divergent Perturbations on a Sphere

Abstract The linear instability of divergent perturbations that evolve on a cos2 mean steady zonal jet embedded in a zonal channel on the β plane and on a rotating sphere is studied for zonally propagating wavelike perturbations of the shallow-water equations. The complex phase speeds result from the imposition of the no-flow boundary conditions at the channel walls on the numerical solutions of the linear differential equations for the wave latitude-dependent amplitude. In addition, the same numerical method is applied to the traditional problem of linear instability of nondivergent perturbations on the β plane where results reaffirm the classical, analytically derived, features. For these nondivergent perturbations, the present study shows that the growth rate increases monotonically with the jet maximal speed and that the classical result of a local maximum at some finite westward-directed speed results from scaling the growth rates on the jet’s speed. In contrast to nondivergent perturbations, divergent perturbations on the β plane have no short-wave cutoff, and so the nondivergent solution does not provide an estimate for the divergent solution, even when the ocean is 1000 km deep (i.e., when the speed of gravity waves exceeds 10 Mach). For realistic values of the ocean depth, the growth rates of divergent perturbations are smaller than those of nondivergent perturbations, but with the increase in the ocean depth they become larger than those of nondivergent perturbations. For both perturbations, a slight asymmetry exists between eastward- and westward-flowing jets. The growth rates of divergent perturbations on a sphere are similar to those on the β plane for the same values of the model parameters, but the asymmetry between eastward and westward jets is more conspicuous on a sphere. The value of g′H′ (g′ is the reduced gravity; H′ is the equivalent mean layer thickness), which is filtered out in nondivergent theory, determines for divergent perturbations the relative magnitude of zonal velocity, meridional velocity, and height but has little effect on the growth rates.

Modeling stream–aquifer interactions with linear response functions

The problem of stream–aquifer interactions is pertinent to conjunctive—use management of water resources and riparian zone hydrology. Closed form solutions expressed as convolution integrals of impulse response and unit step response functions are derived for channel flow and stream–aquifer interactions. The solutions, obtained using Laplace transforms, relate channel reach discharge, stream–aquifer exchange rates, and associated flow volumes to hydrologic processes and regulatory and management control measures. Channel flow is modeled using a simple mass balance equation and the Muskingum linear storage relationship, and groundwater flow is approximated by a linearized representation of the one-dimensional Boussinesq equation. Within a given channel reach, the aquifer is assumed to be homogeneous, separated from the channel by a semipervious layer, and with a time-variable prescribed head or no-flow boundary condition at a finite distance normal to the channel flow direction. Discrete-time unit response functions are derived for arbitrary channel inflow hydrographs and other system's excitations, which extend the utility of the model to a wider spectrum of water management problems. Although the closed-form solutions are based on simplifying assumptions which may limit the utility of the solutions, applications to example flow scenarios nonetheless illustrate the robustness of the solutions to a variety of applications such as the bank storage problem, effect of drought periods, and groundwater–surface water interactions in the presence of water management controls.

Laplace-Domain Solutions for Radial Two-Zone Flow Equations under the Conditions of Constant-Head and Partially Penetrating Well

A mathematical model is presented for a constant-head test performed in a partially penetrating well with a finite-thickness skin. The model uses a no-flow boundary condition for the casing and a constant-head boundary condition for the screen to represent the partially penetrating well. The Laplace-domain solutions for the dimensionless flow rate at the wellbore and the hydraulic heads in the skin and formation zones are derived using the Laplace and finite Fourier cosine transforms. The solutions of hydraulic heads have been shown to satisfy the governing equations, related boundary conditions, and continuity requirements for the pressure head and flow rate at the interface of the skin zone and undisturbed formation. In addition, an efficient algorithm for evaluating those solutions is also presented. The dimensionless flow rates obtained from new solutions have been shown to be better than those of Novakowski's solutions, especially when the penetration ratio is large.

Upscaled models of flow and transport in faulted sandstone: boundary condition effects and explicit fracture modelling

Faults formed by shearing of joint zones in sandstone contain fine-scale features that cannot be represented explicitly in large-scale flow simulations. Upscaled models are, therefore, required for reservoir engineering computations. These models attempt to capture fine-scale effects through equivalent permeabilities that are computed from the underlying fine-scale characterization. In this paper the impact of several different local boundary conditions on the calculated equivalent permeability is assessed. Pressure–no-flow, periodic and mirror-periodic boundary specifications are considered. The resulting coarse-scale permeability tensors are shown to be highly dependent on the local boundary conditions used in the models. In cases with through-going high-permeability features, pressure–no-flow and mirror-periodic boundary conditions provide upscaled permeabilities that correctly capture global flow characteristics. Periodic boundary conditions, by contrast, are more suitable for systems lacking through-going high-permeability features. This sensitivity to boundary conditions calls into question the robustness of the equivalent permeability for the general case and suggests that dominant through-going features would best be modelled explicitly. In addition, due to the very small thickness and high permeability of some through-going structural features (e.g. slip surfaces), globally upscaled models are inadequate for the modelling of transport. To address these issues, a ‘partial upscaling’ method – removing the through-going high-permeability features from the fine model, upscaling to a coarse grid and then reintroducing the high-permeability features back into the coarsened model – is adopted. This procedure is shown to provide coarse models that give accurate predictions for both flow and transport.

Ewald Method for the Analytic Solution of Simple Reservoir Problems With Neumann Boundary Conditions

Summary Analytic solutions to reservoir flow equations are difficult to obtain in all but the simplest of problems. Under the assumption of steady-state conditions, incompressible flow, and a translationally invariant permeability tensor, the Darcy equation leads to Laplace's equation for the hydrostatic potential. Methods of potential theory in both two and three dimensions may then be brought to bear on the problem, which is usually classified by the type of boundary condition. Dirichlet and Neumann-type boundary conditions, where the value of the potential and the flow respectively are specified on the boundary, are elegantly solved by simple superposition of sources in the case where the boundary is at infinity, or by Green's function techniques in cases with finite boundaries. Reservoir simulators, limited by the finite extent of their models, impose conditions on the flow at the boundary, which through Darcy's equation is proportional to the gradient of the potential. Analytic solutions to these Neumann boundary conditions are in general approached using the method of images to explicitly construct the Green's function for a point source. However, as can be easily shown, problems with no-flow boundary conditions lead to an infinite array of images. The sum of the potentials of the infinite lattice results in the solution in the simulation box. In this paper, we make use of the direct analogy between the reservoir flow equations and electrostatics to illustrate a method commonly used in the field of solid-state physics to sum the potentials of crystal structures. The advantage of this method is that it is fast, easily implemented, and may be used for an arbitrary configuration of injecting and producing wells, as long as net flow is zero. In addition, it is useful for a very general class of simulation boxes, not limited to orthogonal axes.

Experimental Study of the Impact of Boundary Conditions on Oil Recovery by Co-Current and Counter-Current Spontaneous Imbibition

Spontaneous imbibition (SI) is a very important oil recovery mechanism from low-permeability fractured reservoirs. The impact of different sample shapes and boundary conditions on oil recovery by counter-current and co-current SI into low-permeability (2 mD) strongly water-wet outcrop chalk have been investigated. The experimental program was divided into two parts. In part 1, counter-current SI was studied by submerging cubic, cylindrical, and irregular rock samples totally in the aqueous phase, whereas co-current SI was investigated in part 2 by exposing different portions of a cubic rock sample to an oil-phase. The scaling law of Ma et al. -1 was used to correlate counter-current SI data with all parameters than the characteristic length, L C , held essentially constant. The results confirm the ability for L C to account for cylindrical (diameters in the range from 2 to 10 cm) and cubic sample shapes and different no-flow surface boundary conditions (2, 4 and 5 faces were coated with polyester on the cubic rock samples; top and bottom faces on the cylindrical cores). Rock samples having irregular shapes (increasing and decreasing cross-sectional areas were investigated) could not be scaled according to L C because the overall shape of the imbibition curves deviated from curves generated by rock samples having regular shape. No scaling law correlating imbibition data for different boundary conditions exists for tests performed under co-current flow conditions; hence the results in part 2 of the experimental program are interpreted qualitatively. When using the SI rate for a cube with all faces open to imbibition as the reference (counter-current flow conditions), the results show that depending on the area fraction of the cube covered by oil, co-current SI can either be faster (water-oil area-ratio (WOAR) > 1), equal (WOAR = 1) or slower (WOAR < 1)) than the reference system. In all cases, ultimate recovery is significantly higher for tests performed under co-current flow conditions than under counter-current (70% of initially oil in place (IOIP) vs 58% IOIP for counter-current). The results indicate that when evaluating oil recovery potential on reservoir rock samples due to SI experimentally, tests should be performed both under co-current and counter-current flow conditions because the latter one might under certain circumstances predict too low rates and ultimate recoveries.

Integrated Geological Modelling and Transmissibility Scale-Up for Improved Accuracy in Reservoir Simulation

Abstract This paper presents a study on the development of efficient methods for scaling petrophysical properties from high-resolution geological models to the resolution of reservoir simulation models. These methods were evaluated using data for the Gypsy field located in northeastern Oklahoma near Lake Keystone. Various criteria that may be used to reduce the number of layers to represent a reservoir in simulation studies were investigated. The approach attempts to reduce the number of grid blocks needed by combining thinner layers and representing them by scaled petrophysical properties assigned to the resulting thicker layers. Three different geological models were developed based on channel identifiers, lithofacies, and flow units, respectively. The effect of the criteria used for combining layers on simulation results was studied by conducting scale-up for three different geological models, three different production scenarios, and three different boundary conditions. In addition to the linear flow scale-up of transmissibility between two grid blocks, a scale-up of the productivity index (PI) was found to be important and necessary in order to account for the radial flow around the wellbore. Special consideration was also needed for the pinch-out grid blocks in the system. The validity of the proposed approach was evaluated by comparing the performance prediction for various reservoir flow scenarios using fine-scale and coarse-scale reservoir models. Strategies of geological modelling were found to have a significant impact on simulation results, especially during the early phase of flow. The use of lithofacies as a criterion provided the closest match to fine-scale results. Amongst various drive mechanisms compared, the best matches for both water production and reservoir pressure were achieved for the line-drive scenario. A better match was obtained for no-flow boundary conditions compared to an open boundary conditions scenario. Introduction Various scale-up techniques have been developed in recent years, such as the averaging method(1, 2), the tensor method(3–6), transmissibility scale-up(7, 8), renormalization technique(9–12), and the pressure-solver method(2, 13). White and Horne(7) showed that the general tensor scaling procedure can give estimates on a coarse grid that are significantly more accurate than other scale-up methods, but greatly increases the computation effort. It is important that the scale-up near wells includes the characteristics of radial flow. Soeriawinata and Kelkar(16) presented an analytical method in which the well block was divided into sectors. The permeability of the well block in the upscaled layers, was determined using a thickness averaging method. Ding(17) calculated the equivalent coarse grid transmissibility based on fine grid simulations and scaled the productivity index using an imposed well condition for radial flow near a well. It was reported that the errors between fine and coarse-scale results, using a scaleup procedure including a radial flow region, are much lower than using a standard procedure. Many scale-up methods concentrate on the mathematics of combining petrophysical properties, but ignore the geological heterogeneity and structural details.

Hydrocarbon Drainage along Corners of Noncircular Capillaries

An approximate solution for liquid flow along corners of noncircular capillaries is proposed that relates the flow resistance to the geometry of a capillary and to the contact angle of the interfaces. The theory predicts liquid flow rates for both two-phase (hydrocarbon/air) and three-phase (hydrocarbon/air/water) drainage. This solution is found to be superior to other expressions in the literature. Drainage rates of two different hydrocarbons along corners of square capillaries are measured with both air and water as stationary phases. The new solution was used to predict the measured hydrocarbon drainage rates successfully. Comparison of the measured and predicted drainage rates indicates that a free boundary is appropriate for the air/hydrocarbon interface and a no-flow boundary condition is valid for the water/hydrocarbon interface. The effect of spreading coefficient on drainage rates is also demonstrated by the measurements and is compared with the proposed solution. A brief discussion of three-phase relative permeabilities is offered based on the measured drainage rates and proposed solution.

Upscaling for Reservoir Simulation

Upscaling for Reservoir Simulation

An Efficient Algorithm For Removal Of Inactive Blocks In Reservoir Simulation

Abstract In the efficient simulation of reservoirs having irregular boundaries one is confronted with two problems: the removal of inactive blocks at the matrix level and the development and application of a variable band-width solver. This paper presents a simple algorithm that provides effective solutions to these two problems. The algorithm is demonstrated for both the natural ordering and D4 ordering schemes. It can be easily incorporated in existing simulators and results in significant savings in CPU and matrix storage requirements. The removal of inactive blocks at the matrix level plays a major role in effecting these savings whereas the application of a variable band-width solver plays an enhancing role only. The value of this algorithm lies in the fact that it takes advantage Of irregular reservoir boundaries that are invariably encountered in almost all practical applications of reservoir simulation. Introduction Reservoir simulation has become a valuable tool and has been effectively used in the development and optimization of oil recovery from petroleum reservoirs. The application of a reservoir simulator to a particular field provides useful information about the reservoir performance under various conditions. To obtain such information, the practicing engineer carries out a process that involves three consecutive tasks - data preparation, history matching, and future performance prediction. The actual reservoir, in this process, is represented in the simulator by a system of equations written for a collection of regular blocks. Almost all practic:a1 applications of reservoir simulation involve reservoirs having irregular boundaries. An example of a two dimensional reservoir is presented in Figure 1. The figure illustrates how the reservoir is subdivided into blocks. Blocks lying on the reservoir boundaries, which are identified by "x", are called "boundary" blocks. Blocks that are within and including boundary blocks are called "active" blocks whereas those exterior to the boundary blocks, which are cross-hatched, are called "inactive" blocks. In practice, all blocks are assigned dimensions and properties such as thickness, elevation, porosity, and permeability in the X-, y-, and z-directions. The dimensions and properties of boundary blocks are modified to ensure a good approximation of the actual reservoir boundaries. The inactive blocks are assigned zero permeability in the three principal directions. This, in effect, represents no-flow boundary condition between the active boundary blocks and neighbouring inactive blocks. Other boundary conditions such as constant flow or constant pressure can then be treated as source (sink) terms in the flow equations for boundary blocks. The simulation of the actual reservoir is achieved by solving all equations for the computational grid. However, we recognize, in this example, two computational grids. The first computational grid (CG1) considers all blocks regardless of their being active or inactive, and as a result it has regular boundaries coinciding with the external boundaries of the mesh itself. The second computational grid (CG2) considers active blocks only, and as a consequence it has irregular boundaries (Fig. 1). Although both computational grids produce identical results, each has its pros and cons.