All examples in the book are based on the MATLAB Reservoir Simulation Toolbox (MRST), an open source toolbox popular in both academic institutions and the petroleum industry. In the direction of the agas cap, the tank contains more gas than can be dissolved in the oil.

## Reservoir Simulation

Upscaling was and probably still is one of the most active areas of research in the oil industry. It is obviously important to continue improving such models to gain a better understanding of the reservoir and reduce uncertainty.

## Outline of the Book

The first part of the book discusses how to represent a geological medium as a discrete computer model that can be used to study the flow of one or more fluid phases. Chapter 2 gives you a crash course in petroleum geology and geological modeling, written by a non-geologist. The fourth and final part of the book is devoted to discussing additional computational methods and tools that can be used to address common tasks within reservoir engineering workflows.

## The First Encounter with MRST

MRST also requires porosity (i.e., the void fraction of the bulk volume) to be given for all models. Unlike grids, petrophysical data and boundary conditions, data structures for representing fluid properties are not part of the basic functionality of MRST.

## Formation of Sedimentary Rocks

A floodplain is the area of land adjacent to the river bank that will be covered by water when a river breaks its bank and floods during periods of high water discharge. Reverse dip-slip faults occur as a result of compressive shortening of the crust, and are sometimes called convergent faults, as rock volumes on opposite sides of the fault plane move horizontally towards each other.

## Creation of Crude Oil and Natural Gas

Stratigraphic traps are particularly important in hydrocarbon exploration; about 4/5 of the world's oil and gas reserves are located in anticlinal traps. The drilling of submerged oil wells began just before the turn of the nineteenth century: from platforms built on piles in the fresh waters of Grand Lake St.

## Multiscale Modeling of Permeable Rocks

Unfortunately, building a geological model for a reservoir is like finishing a puzzle with most of the pieces missing. However, extrapolating properties from the pore scale to an entire reservoir is very challenging, even if the entire pore space of the reservoir were known.

## Modeling Rock Properties

For a completely solid medium, porosity is a static, dimensionless quantity that can be measured in the absence of flow. Here, the diagonal terms {Kxx,Kyy,Kzz} represent how the flow rate in an axial direction depends on the pressure drop in the same direction.

## Property Modeling in MRST

The lower left graph shows the net to gross, i.e., reservoir rock fraction per bulk volume. Although these properties can be calculated from the geometry (and topology) of the network, it is often useful to calculate them in advance and include them explicitly in the network representation.

## Structured Grids

The grid is given by the coordinates of the vertices, but a mapping exists that will transform the curvilinear grid into a uniform Cartesian grid so that each cell can still be referenced using a multi-index(i1,i2, . The active part of the model is shown in yellow color in the right-hand graph. In the figures, the image is taken from penny, one of the standard datasets distributed with MATLAB, and then used to define the geometry of the grid and determine permeability assign values.

## Unstructured Grids

Each vertex in the Delaunay triangulation corresponds to and is the center of a Voronoi cell. To avoid this situation, one can perform a constrained Delaunay triangulation and insert additional points where the constraint is not met (i.e., the center of the perimeter is outside the triangle). Try using a combination of the triangle Grid and Pebit to approximate such a grid, as shown on the right in the image.

## Stratigraphic Grids

In the next two subsections, we discuss the basics of the two most commonly used forms of stratigraphic grids. So far, the topology and geometry of a vertex grid has not deviated from the mapped Cartesian grids studied in the previous section. The original corner point format has been extended in several ways, for example to allow the vertical intersection of two coordinate lines in the form of the letter Y.

## Grid Structure in MRST

One of the entries in neighbors(i,:), but not both, can be zero, to indicate that facei is an external face belonging to only one cell (the non-zero entry). In total, the grid consists of eight cells, fourteen faces and seven nodes, which are shown in Figure 3.29, together with the contents of the fields cells.faces, faces.nodes and faces.neighbors. In the original implementation, all cell quantities were calculated inside a loop, which may not be as efficient as calculating the face quantities.

## Examples of More Complex Grids

Next, let's inspect the error structure in the lower right corner of the plot. As another example, we will generate a Cartesian grid that has a radial refinement around two wells in the interior of the domain. Using this technique, you can more easily see the structure of the model's interior.

## Fundamental Concept: Darcy’s Law

If you have read the chapters of the book in chronological order, you have already encountered the equations that model the flow of a single, irreversible fluid through a porous medium twice: first in section 1.4, which showed how the vertical equilibrium inside can be calculated using MRST, and then in section 2.4.2 we discussed the concept of exchange. We then introduce basic spatial discretizations for the resulting equations and show how these discretizations can be formulated as abstract discrete differentiation operators. A fuller discussion of how these discretizations can be turned into efficient solvers is deferred to Chapters 5 and 6, which focus on the incompressible case, and Chapter 7, which discusses compressible flow and how discrete differentiation operators can be combined with automatic differentiation for rapid prototyping of new solvers.

## General Flow Equations for Single-Phase Flow

To develop the differential equation, we first assume that the porosity and fluid viscosity do not depend on pressure. In the absence of gravitational forces and source terms, this is a linear equation for the fluid density, similar to the classical heat equation with variable coefficients. Here, ω denotes the acentric factor of the species, which is a measure of the centricity (deviation from spherical shape) of the molecules in the liquid.

## Auxiliary Conditions and Equations

The largest percentage of the pressure drop associated with a well occurs near the well and thus the pressure at the well radius will deviate significantly from the volumetric pressure average in the cell. In the above expression, we introduced the formation volume factor B, defined as the ratio between the volume of the fluid at reservoir conditions and the volume of the fluid at surface conditions. For deviated wells, h indicates the length of the grid block in the main direction of the wellbore and not the length of the wellbore.

## Basic Finite-Volume Discretizations

However, in a finite volume method we only know the cell average value of the pressure inside the cell. Let us now construct the discrete analogs of the divergence and gradient operators, which we denote divand degree. In the continuous case, the gradient operator is the adjoint of the divergence operator (up to a sign), which we have.

## Basic Data Structures in a Simulation Model

These two fields are vectors of length equal to the number of completions in the well and contain the bottomhole pressure and flow for each completion. For convenience, fees and may contain a single value; this value is then used for all faces specified in the call. The solvers assume that the boundary conditions are given at the boundary; conditions inside the domain give unpredictable results.

## Incompressible Two-Point Pressure Solver

To focus on the discretization and keep the discussion simple, we will not look at the full implementation of the two-point solver inincomp. We then go through all the planes and calculate the transmissivity of the plane as the harmonic mean of the half transmissivity. In our case, G.cells.faces(:,1) returns the global face number for each half-face, hence the call to accu-array sums the portability of the half-faces corresponding to a given global face and stores the result in the correct place in a vector of G.faces.number elements.

## Upwind Solver for Time-of-Flight and Tracer

As you may recall, time of flight is the time it takes a neutral particle to travel from the nearest fluid source or inlet boundary to any point in the reservoir. You can also calculate the backward time of flight—the time it takes to travel from any point in the reservoir to the nearest liquid sink or outlet boundary—with the same equation if we change the sign of the flow field and modify the boundary conditions and/or source terms accordingly. We have now created the full discretization matrix and the time-of-flight can be calculated by a simple inversion of the matrix.

## Simulation Examples

The left graph shows the pressure distribution and in the right graph we imposed streamlines passing through the centers of the cells on the northwest-southeast diagonal. In particular, we see that the points in the radial plot closely follow those of the reference on a fine scale. The producers are specified in the same way. Figure 5.11 shows the well positions and pressure distribution.

## The TPFA Method Is Not Consistent

In the two-point approximation (6.1) we only use the pressures pi and pk to approximate the flux. To link the derivation in example 6.1.1 to the geometry of the cells, we can use the one-sided definition of fluxes (4.51) of the two-point diagram to write. For the two-point scheme to be correct, the second term in the integrand must be zero.

## The Mixed Finite-Element Method

We can also prove that the solution(p,v) is a saddle point of the Lagrange functional, i.e. that L(v,w) ≤L(v,p) ≤L(u,p) for all. There are four basis functions that have nonzero support in the interior of the center cell. The matrix entries depend on the geometry of the grid cells and whether it is isotropic or anisotropic.