Special topics in upscaling

Quick refresher

The equivalent permeability of a cell is well defined when the heterogeneity is much finer than the cell, and the flow changes over distances much larger than the cell (scale separation). For cells about as large as the correlation length, the answer depends most strongly on the boundary conditions. In the calculation of the lesson, at $\sigma_{\ln k}$ = 1.5 ( $\sigma_{\ln k}$ is the standard deviation of the logarithm of permeability, a measure of how widely it varies), in 56% of the blocks the two kinds of conditions differ by more than 20%.
The most dangerous bodies are thin but long: shale lenses and thin layers, tight layers, conducting thin layers, fractures. Sealed side faces of a cell do not let the flow go round a lens through the neighbouring cell. A uniform pressure on the top and bottom faces acts as a conducting film, and lets the flow go round the lens too easily.
Heterogeneity finer than a cell of the geological model (laminae, cross-bedded sets, which are groups of inclined laminae) is transferred in a separate first step. At the scale of laminae the flow is almost everywhere capillary, and the capillary equilibrium limit is suitable for it. At the scale of a coarse cell the flow regime (which forces dominate) changes from viscous near wells to mixed between them.
The usual two-point scheme sees only the diagonal of the tensor. If the principal axes of permeability (the directions of the largest and the smallest permeability) are turned relative to the grid, the off-diagonal term is lost. With a turn of 45° the cell looks isotropic, whatever the anisotropy is.
Near a well the flow is radial (it converges on the well from all sides). If the rock at the wellbore is heterogeneous, a cell with a well needs a well index (the factor that links the well rate to the pressure drop) found by a local calculation on a fine grid. This index hardly depends on the size of the calculation region, while the ratio of the rates does depend on it. So we correct the index.
The single-phase rate checks only the total. The distribution of the flow over fast and slow paths is shown by the flow capacity – storage capacity curve and by the first arrival of a tracer. More layers do not mean more accuracy: what matters is which layers are combined.

Key relationships

\lambda \ll \Delta \ll L

scale separation: heterogeneity, then cell, then the scale of change of the flow

k(\varphi) = k_{xx}\cos^2\varphi + k_{yy}\sin^2\varphi + 2k_{xy}\sin\varphi\cos\varphi

directional permeability; a scheme without $k_{xy}$ loses the last term

u(r) = \dfrac{q}{2\pi r h}

Darcy velocity near an injection well: $N_{cv}$ and $N_{gv}$ (the ratios of capillary and of gravity forces to viscous forces) grow in proportion to the distance

W_{iD,1} = \dfrac{\Delta C_1}{\Delta F_1}

first arrival of the tracer through the fastest layer, in pore volumes

Common pitfalls

Taking the equivalent permeability of cells as large as the correlation length without comparing it for two kinds of boundary conditions.
Thinking that thinner layers of the coarse grid always give a more accurate vertical flow.
Giving a cell of the geological model the permeability and the relative permeability curves of samples, and skipping the step of laminae.
Dropping the off-diagonal term of the tensor with strong anisotropy and principal axes at 20–70° to the grid.
Correcting the rate of a well by a multiplier chosen during history matching, instead of the index from a local calculation.
Checking a coarsened model only in the conditions for which it was built.

After this lesson you can

explain the scale separation condition, and use the correlation length, the cell size and the distance between wells to estimate how reliable the equivalent permeability of a cell is
show by calculation why discontinuous shale layers and other thin, long bodies are the main source of errors of vertical permeability, and how the sign of the error depends on the shape of the coarse cell
describe stepwise upscaling (upscaling in several steps, from laminae, which are very thin layers, to a cell of the geological model), and choose a method for each step and each part of the reservoir by the balance of forces
explain what the two-point scheme (the usual scheme that finds the flow through a face from the pressures of two cells only) loses with anisotropy at an angle to the grid, and choose between aligning the grid, a multipoint scheme (one that uses more cells around a face) and neglecting the off-diagonal term
calculate the well index from a local calculation on a fine grid, and explain why the correction is put into the index and not into the rate
check a coarsened model by the flow capacity – storage capacity curve (a curve that shows how much of the flow goes through how much of the pore volume) and the first arrival of a tracer (a marker substance carried by the flow), and recognise the cases in which upscaling does not work

Scale separation: when a cell has a permeability of its own

Upscaling replaces a non-uniform block by a uniform one with the equivalent permeability $k^{*}$ (flow calculation in the block). The replacement makes sense under two conditions. First, the cell holds many independent patches of heterogeneity, and its $k^{*}$ depends neither on their random arrangement nor on what lies beyond its faces. Second, the flow inside the cell is almost uniform: the pressure and the saturation change over distances much larger than the cell itself. Together these conditions are called scale separation.

The condition of scale separation

\lambda \;\ll\; \Delta \;\ll\; L

where:

$\lambda$ size of the heterogeneity in the given direction: the correlation length of permeability, m
$\Delta$ size of the cell in the same direction, m
$L$ distance over which the flow changes noticeably: between wells, the width of the displacement front, the zone of radial inflow near a well, m

Units: any equal units of length

Holds when

The condition is checked separately for each direction: the correlation lengths in the horizontal and the vertical directions differ by tens or hundreds of times.
“Much less” means, in practice, a difference of a few times. When $\Delta$ is comparable to $\lambda$ , the equivalent permeability of the cell depends on the boundary conditions of its calculation (the figure below).

The left inequality makes $k^{*}$ a property of the cell. The right one allows us to treat the flow in it as linear, as it is treated in the calculation. The smallest volume in which an averaged property stops depending on the size of the volume itself is called the representative elementary volume (REV). The same idea works at the pore scale (Pore-scale flow modelling).

In a reservoir model both conditions are rarely met. The correlation lengths of permeability are from tens of metres to kilometres horizontally, and from centimetres to metres vertically (correlation length). The cells are 50–200 m across horizontally and 1–10 m thick, which is the same order. Between neighbouring wells there are only a few cells. A calculation on synthetic fields shows how dangerous this is. Blocks of different sizes were cut from random lognormal fields (permeability whose logarithm has a normal distribution: a few high values and many low ones). For each block we found $k^{*}_x$ with sealed side faces and with linear conditions (boundary conditions).

CalculatedBlocks from 25 to 800 m cut from lognormal fields with a correlation length of 150 m (10 random versions, or realisations, of a field of 1.6 × 1.6 km on cells of 25 m). Top: the median difference of

k^{*}_x

between linear conditions and sealed side faces, for three values of

\sigma_{\ln k}

. Bottom: the share of the scatter of

\ln k

that is kept between blocks (at

\sigma_{\ln k}

= 1.5; at 1.0 and 2.0 the curves are almost the same). The difference is largest for blocks of about the correlation length, and it grows quickly with the heterogeneity.

A cell much smaller than the correlation length ( $\Delta/\lambda$ < 0.3). The rock inside the cell is almost uniform. The boundary conditions hardly matter, and the coarse cells keep almost all the scatter of properties. Upscaling changes the model little, and saves little.
A cell of the order of the correlation length ( $\Delta/\lambda$ from 0.5 to 3). The cell holds one or two bodies, which may or may not connect with the neighbouring ones, and $k^{*}$ depends on the assumption about the neighbours. At $\sigma_{\ln k}$ = 1.5, for blocks of 1.3 correlation lengths, the median difference is 23%, and for 56% of the blocks it is more than 20%. This is the least reliable size, and the most common one.
A cell much larger than the correlation length ( $\Delta/\lambda$ > 5). The heterogeneity inside the cell averages out, $k^{*}$ settles, and the boundary conditions matter less and less. But the scatter between cells also almost disappears: at 5.3 correlation lengths, 9% of the original scatter of $\ln k$ remains between the blocks. The heterogeneity of the scale $\lambda$ has gone from the model into the cells. For single-phase flow this is right. For displacement it is not: water still goes along the best bodies inside the cell, and we put this back with pseudo functions (Two-phase upscaling: pseudo-functions).

The difference grows quickly with the heterogeneity. For blocks of 1.3 correlation lengths the median is 9% at $\sigma_{\ln k}$ = 1.0, 23% at 1.5 and 48% at 2.0. The calculation is for an isotropic plane field. Elongated bodies and three-dimensional flow change the numbers but not the picture (see the advanced section).

Worked exampleIs there scale separation in a model of a river-deposit reservoir

Given (assumed data)

Correlation length of permeability: along the channel belts 1200 m, across them 300 m, vertically 1.5 m
Heterogeneity: $\sigma_{\ln k}$ ≈ 1.5
Cells of the model: 100 × 100 m, layers of 2 m
Distance between wells: 500 m

Solution

The ratio of the cell size to the correlation length: along the belts $100/1200 = 0.083$ , across them $100/300 = 0.33$ , and vertically $2/1.5 = 1.33$ .
From the figure, at $\sigma_{\ln k}$ = 1.5: at the ratio 0.33 the median difference between the boundary conditions is 5.5%, and it is more than 20% for 21% of the blocks. At the ratio 1.33 it is 23%, and more than 20% for 56% of the blocks. Along the belts the difference is even smaller than across them.
Between wells there are $500/100 = 5$ cells. Near a wellbore the pressure changes noticeably already from one cell to the next, and the right-hand inequality is violated.
Layers of 0.5 m would give the vertical ratio $0.5/1.5 = 0.33$ , a difference of 5.5%, but there would be four times more layers.

Answer: Horizontally, the cells are smaller than the correlation length. The grid itself represents the heterogeneity of this scale, and the median difference between the boundary conditions is small. Vertically, the layers of 2 m are comparable with the correlation length. There is no scale separation, and the vertical permeability of the cells is the least reliable property of the model. Near wells the right-hand inequality is also violated.

What it means. The cell size is usually chosen for the horizontal plane and not checked vertically. Here we need either thinner layers where the vertical flow matters (near bottom water and a gas cap, and with gravity separation of the phases), or a check of the vertical transmissibility (how easily fluid flows up and down between two cells) by a calculation on a fine grid (next section). Near wells, a well index of its own is needed (the section on wells). The curves of the figure are for an isotropic plane field and serve as a guide, not as a norm. The data in this example are assumed for illustration.

Thin long bodies: lenses, thin layers, fractures

The correlation length and the histogram describe the heterogeneity, but they do not say the main thing for the flow: do the good patches join into continuous paths, and do the bad ones block the path of the flow? Models with the same histogram and variogram can be connected in quite different ways. Connectivity is set by the shape and the extent of the bodies (connectivity and percolation). The two-point statistics, the variogram (a curve that shows how fast the similarity of two values falls as the distance between them grows), do not determine it.

For upscaling the most dangerous bodies are thin but long: shale lenses and thin layers, tight carbonate layers, thin layers of coarse sandstone, open fractures, fault zones. They are thinner than a cell, so the cell averages them. And they are longer than a cell, so their action is not local: a flow stopped by a lens in one cell goes round it through the neighbouring cell. We take the most frequent case, vertical flow through sandstone with discontinuous shale lenses (baffles, bodies that obstruct the flow).

CalculatedOne random version of the cross-section of the example: sandstone of 100 mD, 1000 m long and 10 m thick (the vertical scale is exaggerated 25 times), with shale lenses 0.5 m thick and about 100 m long, making 10% of the volume. The vertical lines are the boundaries of coarse cells 100 m wide. Many lenses cross them, and the flow goes round a lens through the neighbouring cell.

Worked exampleShale lenses: the vertical permeability on the fine and the coarse grids

Given (assumed data)

Cross-section: 1000 × 10 m; fine grid of 200 × 20 cells of 5 × 0.5 m; sandstone of 100 mD
Lenses: thickness 0.5 m, permeability 0.001 mD, about 10% of the volume; length $\ell$ = 10, 100 or 400 m, random position; 12 realisations for each length
Coarse grids: columns of 100 m; one layer of 10 m, two layers of 5 m, or twenty of 0.5 m; in each cell, $k^{*}_z$ and $k^{*}_x$ come from a flow calculation with sealed side faces
Check: vertical flow through the whole cross-section: the base and the top at different constant pressures, the ends closed

Solution

Fine grid, the average $k_v$ of the cross-section: 24.2 mD at $\ell$ = 10 m, 10.6 mD at 100 m, and 7.0 mD at 400 m. The rule “share of sandstone × its permeability” would give 90 mD for any length.
Coarse grid of one layer, the median of the ratio of its $k_v$ to the fine one: 0.98, 0.98, 0.93. At $\ell$ = 400 m the single realisations range from 0.08 to 1.00.
Two layers of 5 m: 1.53, 1.51, 1.07. At $\ell$ = 400 m the range is from 0.14 to 11.7.
Twenty layers of 0.5 m, that is, coarsening horizontally only: 3.7, 5.0, 1.9.

Answer: With almost the same share of shale, the

k_v

of the cross-section falls from 24.2 to 7.0 mD as the lenses get longer. The coarse grid of one layer is close to the fine one in the median, but when the lenses are longer than the cell it is wrong by several times in some realisations. The grids with inner faces overestimate

k_v

, by up to 5.0 times at

\ell

= 100 m.

What it means. The two opposite errors are set by the way the cell calculation itself is posed. Sealed side faces do not let the flow go round a lens through the neighbouring column. This underestimates $k_v$ , and most strongly when a lens covers the column completely, which often happens with long lenses. A constant pressure on the top and bottom faces of a cell is equal to a thin conducting film over its whole width. A flow that has passed one gap between lenses spreads freely to any gap above it. This overestimates $k_v$ , in the same way as conducting partitions give the upper Cardwell–Parsons bound (bounds). In one layer of 10 m there are no inner faces, and only the first error remains. In twenty layers there is a face above every row, and the result approaches the rule “share of sandstone × permeability”. Which error wins depends on the length and the share of lenses, on the shape of the cell, and on the random placement of the lenses. The data in this example are assumed for illustration.

CalculatedThe ratio of the vertical permeability of the coarse grid to that of the fine one, against the length of the lenses divided by the width of a coarse cell. The points are medians over 12 realisations, and the bars run from the smallest to the largest value. While the lenses are shorter than the cell, the one-layer grid is almost exact. The grids with inner faces overestimate

k_v

at any length. When the lenses are longer than the cell, the scatter between realisations reaches orders of magnitude.

What is done about thin long bodies

We check the vertical transmissibility of units with lenses and discontinuous thin layers by a flow calculation on a fine grid, for representative cross-sections. If necessary, we choose multipliers for the vertical transmissibility of the coarse model.
Continuous barriers longer than a cell are set explicitly: by zero transmissibility of the faces, or by a separate thin layer. We do not average them in the cell (transmissibility upscaling).
Conducting thin layers are kept as separate layers of the coarse grid (non-uniform grid).
A network of open fractures is a separate medium. It is described by a dual porosity model (the fractures conduct, and the matrix stores the oil; Fractures and naturally fractured reservoirs), and not by the equivalent permeability of the matrix.

Finer than a cell of the geological model: stepwise upscaling

A cell of the geological model, tens of metres across and a fraction of a metre thick, is also non-uniform. It holds laminae, cross-bedded sets, shale drapes (thin films of clay) and thin shale layers a centimetre thick (Primary sedimentary structures). The data about the rock come from an even smaller scale: a core plug is a few centimetres, and a probe permeameter (a small device that measures permeability at a point) on a cut face of the core reads millimetres. So before we transfer properties from the geological grid to the simulation grid, we must transfer them from the sample to a cell of the geological model. This is the first step of upscaling.

This step is often skipped, and the cell gets the permeability and the relative permeability curves of samples, as if it were uniform. There are two mistakes. For one phase, plugs are drilled from the strong laminae (the ones that survive drilling). So the vertical permeability of the cell comes out too high, sometimes by orders of magnitude (a missed thin layer). For two phases, capillary forces in a layered cell redistribute the water among the laminae. So the pseudo relative permeabilities along and across the laminae differ by several times (laminae at capillary equilibrium). Samples do not show this.

So we transfer the properties in several steps, by geologically meaningful scales. The approach is called stepwise upscaling, or the method of geopseudo functions (geopseudos). We build models of typical structures: ripples (small waves of sand), a cross-bedded set, a unit with thin layers. We give the laminae a permeability from point measurements, and relative permeability and capillary pressure curves from core experiments (Petrophysical input for simulation: relative permeability and capillary pressure). We calculate the flow, and we get the equivalent permeability and the pseudo functions for the next step.

Steps of upscaling

Step: Laminae → a layer or a set
Scale: mm → dm
What controls the flow: capillary forces
How it is calculated: the capillary equilibrium limit

Step: Layer → a cell of the geological model
Scale: dm → tens of m × fractions of m
What controls the flow: capillary and viscous forces
How it is calculated: flow calculation in a block, pseudo functions

Step: Cell of the geological model → simulation cell
Scale: → hundreds of m × metres
What controls the flow: viscous forces, gravity, connectivity of bodies
How it is calculated: single-phase upscaling, non-uniform grid, dynamic pseudo functions

Step	Scale	What controls the flow	How it is calculated
Laminae → a layer or a set	mm → dm	capillary forces	the capillary equilibrium limit
Layer → a cell of the geological model	dm → tens of m × fractions of m	capillary and viscous forces	flow calculation in a block, pseudo functions
Cell of the geological model → simulation cell	→ hundreds of m × metres	viscous forces, gravity, connectivity of bodies	single-phase upscaling, non-uniform grid, dynamic pseudo functions

The boundaries of the steps and the regimes are a guide. We check the regime of each step by the balance of forces (the example below, and the section on the balance of forces).

Worked exampleLaminae at capillary equilibrium: up to what velocity

Given (assumed data)

Laminae: from the example of the lesson Two-phase upscaling: pseudo-functions: 400, 60 and 6 mD, thickness about 1 cm; the vertical permeability of the unit is 24.8 mD
Length over which the saturation equalises: 0.3 m, a cell of the detailed model
Fluids: oil 3.0 mPa·s, $k_{ro}'$ = 0.85; capillary pressure difference between the laminae 0.04 MPa
Injection well: 300 m³/d at reservoir conditions, reservoir thickness 16 m

Solution

The ratio of capillary to viscous forces $N_{cv}$ (formula) is equal to one at a Darcy velocity $u^{*} = 0.0853\cdot24.8\cdot0.85\cdot0.04\cdot0.3/(3.0\cdot0.01^2) = 71.9$ m/d.
Near an injection well, $u = q/(2\pi rh)$ is larger than $u^{*}$ only at $r < r^{*} = 300/(2\pi\cdot16\cdot71.9) = 0.042$ m, that is, inside the wellbore.
For laminae 5 cm thick, $u^{*}$ is 25 times smaller: 2.88 m/d, and $r^{*} = 300/(2\pi\cdot16\cdot2.88) = 1.04$ m.

Answer: Laminae a centimetre thick stay at capillary equilibrium everywhere in the reservoir. Laminae 5 cm thick leave it only within 1.04 m of the well.

What it means. At the first step, pseudo functions can be calculated in the capillary equilibrium limit, and they hardly depend on the injection rate. As the layers get thicker, the picture changes quickly. The time of capillary equalisation grows as the square of the thickness, and for units metres thick, capillary equilibrium is no longer reached at a field rate (laminae and the cross-section). There we need dynamic methods. The estimate gives an order of magnitude and not a boundary. The data in this example are assumed for illustration.

The first step is not needed everywhere. It is needed when:

immiscible phases (phases that do not mix) are displaced, by water or gas (for single-phase flow, averaging the permeability with a check of the vertical one is usually enough)
the laminae differ in permeability by several times or more, and so also in capillary pressure
the flow goes across the laminae: vertically, to bottom water and a gas cap, with gravity separation of the phases, or across inclined cross-bedded sets (cross-bedded set)
the fine layering is continuous over the whole reservoir, and its influence is not hidden by larger errors (if the reservoir is a set of disconnected bodies, we first deal with the connectivity of the bodies, and only then look at the small scale in more detail: types of architecture)

Anisotropy at an angle to the grid: what the two-point scheme loses

If the principal directions of permeability (the directions of the largest and the smallest permeability) do not coincide with the grid axes, the equivalent permeability of a cell is a full tensor. This happens with the inclined laminae of cross-bedded sets (derivation for inclined laminae), with beds that dip relative to a horizontal Cartesian grid (a grid of rectangular cells in straight rows and columns), and with channels that run at an angle to the grid lines in the horizontal plane. The block of the lesson on permeability upscaling is such a case: there the periodic conditions gave an off-diagonal term of 20.0 mD. The question of this section is what a simulator does with such a tensor.

The usual scheme of a simulator calculates the flow through a face from the pressures of two cells, with one transmissibility of the face (cell equation). This is the two-point flux approximation (TPFA). It is exact only if the grid is K-orthogonal: the vector $\mathbf{K}\mathbf{n}$ , where $\mathbf{n}$ is the normal to the face (the direction perpendicular to it), points along the line that joins the centres of the two cells. In plain words, the flow through the face then depends only on the pressures of those two cells. A Cartesian grid with a diagonal tensor is K-orthogonal. With a full tensor it is not. A corner-point grid (a grid whose cells are bent to follow the layers and faults) with skewed cells is not either, even with a diagonal tensor (corner-point geometry). Then the two-point scheme loses the flow that is driven by the component of the gradient along the face, and this error does not disappear when the cells are made smaller (see the advanced section).

Directional permeability of a tensor

k(\varphi) = \mathbf{n}^{\mathsf{T}}\mathbf{K}\,\mathbf{n} = k_{xx}\cos^2\varphi + k_{yy}\sin^2\varphi + 2k_{xy}\sin\varphi\cos\varphi

where:

$\varphi$ direction: the angle from the $x$ axis to the $y$ axis, °
$\mathbf{n}$ unit vector of this direction, $(\cos\varphi, \sin\varphi)$ , —
$k_{xx},\ k_{yy},\ k_{xy}$ components of the equivalent permeability tensor in the axes of the grid, mD
$k(\varphi)$ component of the Darcy velocity along $\mathbf{n}$ at a unit pressure gradient along $\mathbf{n}$ , written as a permeability, mD

Units: permeabilities in any equal units

Holds when

The flow is in the $x$ – $y$ plane, and the tensor is symmetric.
The largest and smallest values, the principal permeabilities $k_1$ and $k_2$ , are reached along the principal axes, turned by the angle $\alpha$ : $\tan2\alpha = 2k_{xy}/(k_{xx} - k_{yy})$ .
A scheme that leaves out $k_{xy}$ sees $k_{xx}\cos^2\varphi + k_{yy}\sin^2\varphi$ . Along a principal axis it is wrong by $-(1 - k_2/k_1)\sin^2 2\alpha/2$ .

At $k_{xy}$ = 0 the principal axes coincide with the grid axes, and there is nothing to drop. With principal axes at 45° to the grid, $k_{xx} = k_{yy} = (k_1 + k_2)/2$ . Without the off-diagonal term, the cell looks isotropic, however strong the anisotropy is.

Worked exampleA block with a channel at an angle: what the grid sees without the off-diagonal term

Given (assumed data)

Tensor of the block (periodic conditions): $k_{xx}$ = 62.4, $k_{yy}$ = 64.2, $k_{xy}$ = 20.0 mD, from the lesson Permeability upscaling for single-phase flow

Solution

Principal axes: $\tan2\alpha = 2\cdot20.0/(62.4 - 64.2) = -22.22$ . The numerator is positive and the denominator is negative, so $2\alpha$ is in the second quadrant: $2\alpha = 92.6°$ and $\alpha = 46.3°$ , along the channel.
Principal permeabilities: $(62.4 + 64.2)/2 \pm \sqrt{0.9^2 + 20.0^2} = 63.3 \pm 20.0$ . That is, $k_1$ = 83.3 mD along the channel and $k_2$ = 43.3 mD across it.
The grid without $k_{xy}$ sees, along the channel, $62.4\cdot0.4773 + 64.2\cdot0.5227 = 63.3$ mD. Across the channel it sees $62.4\cdot0.5227 + 64.2\cdot0.4773 = 63.3$ mD. Here $\cos^2\alpha = 0.4773$ and $\sin^2\alpha = 0.5227$ .
The error along the channel is $63.3/83.3 - 1 = -24$ %, and across it is $63.3/43.3 - 1 = +46$ %. The anisotropy $k_1/k_2$ = 1.92 turns into 1.00.

Answer: Without the off-diagonal term the grid sees the block as almost isotropic, about 63.3 mD in all directions, instead of 83.3 mD along the channel and 43.3 mD across it.

What it means. The channel runs almost exactly along the diagonal. This is the worst case for the two-point scheme: $k_{xx}$ and $k_{yy}$ are almost equal, and all the anisotropy sits in $k_{xy}$ . If such bodies cross the model at an angle over hundreds of cells, the two-point model will underestimate the conductivity along them by a quarter, and overestimate it across them by almost a half. Water in the model will go wider and slower than in the reservoir. The data in this example are assumed for illustration.

CalculatedDirectional permeability of the block of the example: the full tensor (solid line) and the same cell without the off-diagonal term (dashed line). The full tensor gives 83.3 mD along the channel and 43.3 mD across it. Without

k_{xy}

the curve is almost flat.

What is done about this:

Align the grid with the structure. The axes of a Cartesian grid are set along the channel belts or along the dip of cross-bedded sets, and a corner-point grid is built along the layers. In the axes of the layer the tensor is almost diagonal. The price is skewed (slanted) cells of the corner-point grid, where the two-point formula also loses accuracy. So we avoid strong skew.
Use a multipoint scheme. The multipoint flux approximation (MPFA) finds the flow through a face from the pressures of several cells. In a plane problem the cell equation links the cell with eight neighbours, and in a three-dimensional one with up to twenty-six. It keeps the full tensor. The price is more coefficients in the equations and a slower solution. With very strong anisotropy or very skewed cells, the pressure may oscillate in a way that is not real (unphysical oscillations).
Check that the loss is small. By the formula above, the error along a principal axis at a turn of 10° and a tenfold anisotropy is about 5%, and at 45° and only a twofold anisotropy it is 25%. A simple check is to repeat the calculation on a grid rotated by 45° and compare (grid orientation effect).

Flow direction and wells

The equivalent permeability is calculated for linear flow along an axis of the cell, and under an assumption about the neighbours. In a reservoir the flow is not like that. Near wells it is radial. Between an injection well and a production well it spreads out and then converges. In an areal pattern element it is diagonal to the grid. When wells are converted to injection and new wells are drilled, the flow changes direction. Upscaling done for one picture of the flow can be wrong for another.

The methods differ in where they take the boundary conditions from. Local methods calculate each cell separately, with sealed, periodic or linear conditions, or with a buffer zone (boundary conditions). Global methods take the conditions from a single-phase calculation of the whole detailed model with wells. The coarse model then “knows” the flow paths for exactly this pattern of wells (global boundary conditions). Local–global methods alternate the two. The calculation of the coarse model sets the conditions at the boundaries of the local calculations, and the local calculations correct the transmissibilities. This goes on until the results stop changing. The closer the conditions are to the real flow, the more accurate the result. But the result is also more strongly tied to the pattern of wells for which it was found.

A cell with a well is a special case. The flow converges on the wellbore, the resistance of each zone is proportional to the logarithm of the ratio of its radii, and the rock at the wellbore matters far more than its share of the volume. An estimate from the skin factor of the zone around the wellbore is in the lesson on permeability upscaling. It works as long as the changed zone is round and lies inside the equivalent radius of the cell. Otherwise the well index is found by near-well upscaling, a local calculation on a fine grid.

Upscaling near a well

1
Cut out the region
Take the cell with the well and a ring of neighbouring cells, on the fine grid.
2
Solve the inflow on the fine grid
Solve single-phase steady flow. The well sits in its own fine cell, with the Peaceman index (the standard formula for the well index of a cell) for the permeability of that cell (Well model in a grid block). On the outer boundary set a constant pressure, or the pressure from a global calculation.
3
Solve the same problem on the coarse grid
Use the upscaled permeabilities and the same condition at the boundary. The well index is not known yet.
4
Choose the index
Find $\mathrm{WI}^{*}$ at which the coarse model gives the same rate at the same pressures. If needed, correct in the same way the transmissibilities between the cell with the well and its neighbours.
5
Check
Repeat with a larger region. The index should hardly change.

Worked exampleA well in a channel: the index from a local calculation

Given (assumed data)

Cells: 100 × 100 m; a vertical well in the centre of the cell, $r_w$ = 0.1 m, wellbore skin 0; the calculation is per 1 m of thickness
Rock: a channel 20 m wide with a permeability of 300 mD passes through the well from west to east and continues in the neighbouring cells; around it the rock is 30 mD
Local region: 3 × 3 coarse cells (300 × 300 m); the fine grid of 4-m cells has 5625 cells; the outer boundary is at constant pressure

Solution

The permeability of the coarse cell with the well is exact for a strip. Along the channel $k^{*}_x = (20\cdot300 + 80\cdot30)/100 = 84.0$ mD, and across it $k^{*}_y = 100/(20/300 + 80/30) = 36.6$ mD.
Peaceman index for an anisotropic cell (formula): $r_0$ = 20.2 m, and $\mathrm{WI} = 2\pi\sqrt{84.0\cdot36.6}/\ln(20.2/0.1) = 2\pi\cdot55.4/5.308 = 65.6$ mD·m.
At the same pressure difference between the outer boundary and the bottom hole, the fine grid gives a rate 2.40 times larger than the coarse one.
The index at which the coarse model gives the rate of the fine one is $\mathrm{WI}^{*}$ = 345 mD·m. This is 5.3 times more than the Peaceman index.
Control. For a uniform rock of 30 mD, the local calculation gives 35.3 mD·m against 35.6 mD·m by Peaceman's formula. For a region of 5 × 5 cells with the channel, $\mathrm{WI}^{*}$ = 351 mD·m, and the ratio of the rates is 2.12 instead of 2.40.

Answer: The well index in the channel must be increased 5.3 times: 345 mD·m instead of 65.6 mD·m per metre of thickness.

What it means. The well is in rock of 300 mD, but the coarse cell knows about it only through the averages 84.0 and 36.6 mD, calculated for linear flow. The correction goes into the index and not into the rate. The index is a property of the neighbourhood of the well and hardly depends on the size of the region (345 and 351 mD·m). The ratio of the rates does depend on it, because it includes the resistance of the distant cells, which the coarse model passes on correctly. An estimate from the skin of a round zone does not suit this case: the strip is not round, and it goes far beyond the equivalent radius of 20.2 m. The uniform control shows that the method itself is wrong by about 1%. The data in this example are assumed for illustration.

SchematicThe local region of the example: nine coarse cells of 100 m, a 300 mD channel 20 m wide through the well in the centre, and 30 mD around. The coarse grid sees the channel only as the anisotropy of the cells of the middle row (84.0 and 36.6 mD). The local calculation gives the well an index of 345 mD·m instead of 65.6 mD·m.

A layered cell, with the well opening all the layers. The index of a coarse layer is proportional to the sum $k_ih_i$ of its fine layers. This is the same arithmetic mean as for linear flow along the layers, and permeability upscaling is right here (layers with radial inflow).
Part of a coarse layer is opened. The index is calculated from $k_ih_i$ of the opened fine layers, and not from the share of the thickness. A perforation in the best layer, with the average for the coarse layer, underestimates the index. In the worst layer it overestimates it. The inflow from the layers that are not opened adds a skin of partial penetration (extra resistance, because the well opens only a part of the thickness; Partial penetration and well deviation).
A horizontal wellbore. Its index depends on the permeabilities in the plane perpendicular to the wellbore (index of a horizontal well), and on the fine layer in which the wellbore lies. In a coarse cell this knowledge is lost.
Two phases near a well. Cones of water and gas (the rise of water or gas towards the well in a cone) develop within metres of the wellbore (Water and gas coning). On a coarse grid they are represented by local grid refinement or by a radial grid near the well (local grid refinement). More rarely they are represented by pseudo functions calculated for the well from a cone calculation on a fine radial grid.

The balance of forces: from laminae to the cell, from the well to the middle between wells

Pseudo functions depend on where the water is in the cell, and this is set by the balance of viscous, capillary and gravity forces (steady-state methods). Two ratios of forces are already familiar: capillary to viscous, $N_{cv}$ (formula), and gravity to viscous, $N_{gv}$ (formula). Both are inversely proportional to the Darcy velocity and grow with the cell length. $N_{cv}$ is also inversely proportional to the square of the thickness over which the saturation equalises. So the regime changes both with the scale and with the distance from the well.

Darcy velocity near an injection well

u(r) = \frac{q}{2\pi r h}

where:

$u$ Darcy velocity (the rate per unit area) at the distance $r$ from the axis of the well, m/d
$q$ injection rate of the well at reservoir conditions, m³/d
$r$ distance from the axis of the well, m
$h$ thickness of the reservoir that takes the water, m

Units: m, d, m³ (any consistent units)

Holds when

Plane radial flow in a unit of constant thickness. The distribution of the inflow among the layers is not taken into account.
Near other wells the picture of the flow is more complex. For an estimate of the order of magnitude, the formula is applied up to the middle of the distance between wells.

The velocity at the wellbore is hundreds of times larger than in the middle between wells, and $N_{cv}$ and $N_{gv}$ change in the same proportion. A regime that is right for one part of the reservoir is wrong for another.

Worked exampleThe balance of forces along the path from an injection well

Given (assumed data)

Well: 300 m³/d at reservoir conditions, thickness 16 m, 500 m to the production wells
Scale of laminae: the previous example: $N_{cv} = u^{*}/u = 71.9/u$ , with $u$ in m/d
Scale of a coarse cell: length 100 m, thickness 16 m, $k_v$ = 25.8 mD, the cross-section of the lesson Displacement in heterogeneous reservoirs without thin shale layers; oil 3.0 mPa·s, $k_{ro}'$ = 0.85, $\Delta\rho$ = 0.24 g/cm³; capillary pressure difference 0.04 MPa

Solution

Darcy velocity: $u = 300/(2\pi r\cdot16) = 2.984/r$ m/d. This is 2.984 at a distance of 1 m, 0.2984 at 10 m, 0.02984 at 100 m, and 0.01194 at 250 m, in the middle between wells.
Laminae: $N_{cv}$ = 24.1, 241, 2410, 6022.
Coarse cell: $N_{gv} = 0.0853\cdot25.8\cdot0.85\cdot0.00981\cdot0.24\cdot100/(3.0\cdot16\cdot u) = 0.009175/u$ . The values are 0.0031, 0.031, 0.31, 0.77.
For the same cell $N_{cv}/N_{gv} = 0.04/(0.00981\cdot0.24\cdot16) = 1.06$ : the capillary forces are of the same order as gravity.

Answer: The laminae are at capillary equilibrium along the whole path (

N_{cv}

from 24.1 to 6022). The coarse cell is in the viscous regime near the well (

N_{gv}

≈ 0.031 at 10 m), and in the middle between wells all three forces are comparable (

N_{gv}

≈ 0.77).

What it means. Pseudo functions of coarse cells calculated at one rate are right only in part of the reservoir. Near wells the water goes ahead along the permeable layers, as in a viscous regime. Between wells it has time to sink to the base and to redistribute itself by capillary forces. So the tables of pseudo functions are split by zone, near wells and between them. Or we keep the resolution in thickness, so that gravity and capillary forces are calculated on the grid. The data in this example are assumed for illustration.

CalculatedRatios of forces along the path from the injection well of the example (both axes are logarithmic): capillary to viscous at the scale of laminae, and gravity to viscous at the scale of a coarse cell. The capillary forces of the cell almost coincide with gravity and are not shown. The horizontal line shows equal forces.

Which pseudo functions where: steady-state or dynamic

Regime in the cell: capillary
Sign: $N_{cv}$ ≫ 1
Method: the capillary equilibrium limit

Regime in the cell: gravity
Sign: $N_{gv}$ ≫ 1, the layers communicate
Method: vertical equilibrium

Regime in the cell: viscous, the cell is uniform
Sign: $N_{cv}$ and $N_{gv}$ ≪ 1
Method: core curves with a correction for numerical dispersion

Regime in the cell: viscous, the cell is layered or a mosaic
Sign: the same
Method: dynamic pseudo functions or a non-uniform grid

Regime in the cell: intermediate
Sign: ratios of about one
Method: dynamic ones with gravity and capillary pressure, or a finer grid

Regime in the cell	Sign	Method
capillary	$N_{cv}$ ≫ 1	the capillary equilibrium limit
gravity	$N_{gv}$ ≫ 1, the layers communicate	vertical equilibrium
viscous, the cell is uniform	$N_{cv}$ and $N_{gv}$ ≪ 1	core curves with a correction for numerical dispersion
viscous, the cell is layered or a mosaic	the same	dynamic pseudo functions or a non-uniform grid
intermediate	ratios of about one	dynamic ones with gravity and capillary pressure, or a finer grid

The table is a guide to choosing the method, not a calculated result. Steady-state methods are cheap and stable, but right only in their own limit. Dynamic ones cover any regime, but only the one for which they were calculated (where pseudo functions stop working).

Checking: the distribution of the flow, not only its sum

Checks of single-phase flow, such as the rate through a sector and the productivity of wells (how to check a coarsened model), are necessary. But they look at the sum. Displacement depends on how the flow is distributed: what share of it goes along the fast paths, and what share of the pore volume these paths occupy. This picture is given by the flow capacity – storage capacity curve (F–Φ). We order the parts of the reservoir from the fastest to the slowest, and plot the cumulative share of the flow against the cumulative share of the pore volume. For a layered cross-section in single-phase flow this is the Lorenz curve of the layers, a curve that shows how unevenly the flow is shared among the layers (measures of heterogeneity). In the general case we build it from the travel time (see the advanced section).

First arrival of a tracer in a layered cross-section

W_{iD,1} = \frac{\Delta C_1}{\Delta F_1} = \frac{\phi_1h_1\big/\sum_j\phi_jh_j}{k_1h_1\big/\sum_j k_jh_j}

where:

$W_{iD,1}$ the volume injected, in pore volumes of the cross-section, when the tracer arrives through the fastest layer, —
$\Delta C_1,\ \Delta F_1$ shares of the storage capacity and of the flow capacity of this layer: the first point of the F–Φ curve, —
$\phi_j,\ k_j,\ h_j$ porosity, permeability and thickness of layer $j$ ; layer 1 has the largest $k/\phi$ , fraction, mD, m

Units: dimensionless; $k$ and $h$ in any equal units

Holds when

The layers are continuous between the rows of wells, the flow in each is linear, and the pressure difference is the same. In single-phase flow there is then no crossflow between the layers.
The tracer moves with the water without lag and smearing. The rate of each layer is proportional to $k_jh_j$ , and its pore volume to $\phi_jh_j$ .

The steeper the beginning of the F–Φ curve, the earlier the water arrives. A uniform cross-section gives $W_{iD,1}$ = 1: the tracer arrives after one pore volume has been injected.

Worked exampleWhich layers to combine: a check by the F–Φ curve

Given (assumed data)

Cross-section: the eight layers of the lesson Reservoir heterogeneity and permeability averaging: 16 m, from 8 to 900 mD; the sum of $\phi h$ is 2.95 m, and the sum of $kh$ is 2234 mD·m
Coarsening A: four layers, by pairs of layers: 1–2, 3–4, 5–6, 7–8
Coarsening B: three layers: 1–3, the 900 mD layer on its own, 5–8; the non-uniform grid of the lesson Two-phase upscaling: pseudo-functions
Properties of the coarse layers: $k$ and $\phi$ are the averages weighted by thickness

Solution

Coarsening A, the layers (thickness, $k$ , $\phi$ ): 3.5 m, 158.6 mD, 0.2014; 4.0 m, 315 mD, 0.22; 4.5 m, 43.3 mD, 0.1678; 4.0 m, 56 mD, 0.1525.
In decreasing order of $k/\phi$ they are layers 2, 1, 4, 3. The cumulative shares of storage capacity $C$ are 0.298, 0.537, 0.744, 1, and the cumulative shares of flow capacity $F$ are 0.564, 0.812, 0.913, 1.
The Lorenz coefficient (one number for how uneven the flow is) is twice the area between the broken line and the diagonal: 0.34. First arrival: $0.298/0.564 = 0.53$ pore volumes.
Coarsening B: $C$ = 0.085, 0.537, 1 and $F$ = 0.403, 0.812, 1. The Lorenz coefficient is 0.42, and the first arrival is $0.085/0.403 = 0.21$ .
The fine model of eight layers has a Lorenz coefficient of 0.56 and the same first arrival as B. One layer gives 0 and 1.

Answer: The three layers of B kept 75% of the heterogeneity by the Lorenz coefficient, and the first arrival of the fine model. The four layers of A kept 61%, and their tracer arrives 2.5 times later.

What it means. More layers do not mean more accuracy: what matters is which layers are combined. In coarsening A the 900 mD layer merged with its neighbour of 120 mD, and the fast path disappeared. The check takes seconds, and it is done before the displacement calculations. Its conclusion is confirmed by the waterflood calculation in the lesson Two-phase upscaling: pseudo-functions, where grid B reduced the overestimate of production after one pore volume was injected from 16 to 6%. The single-phase F–Φ curve does not see gravity and capillary forces. We check those by a displacement calculation. The data in this example are assumed for illustration.

CalculatedFlow capacity – storage capacity curves of the cross-section of the example: the eight layers, coarsening B (three layers) and A (four layers), and one layer, the diagonal. Curve B repeats the beginning of the fine model, and curve A cuts it off.

What else to check, besides the volumes and the single-phase rate:

the sensitivity of $k^{*}$ to the boundary conditions, a sign of cells as large as the correlation length (the first section)
the vertical communication, by calculating the vertical flow on the fine and the coarse grid for representative cross-sections (the section on lenses)
the loss of off-diagonal terms, by a calculation on a grid rotated by 45°
the well indexes, by a local calculation near typical wells
the displacement, by calculating a sector on the fine and the coarse grid at the field rate, and then at another rate and with another pattern of wells (checking pseudo functions)
the convergence with the grid, by comparison with a model whose cells are twice as small (how many cells are needed)

When there are many models, as in an uncertainty study, we count the upscaling error as one more source of uncertainty. We run a part of the cases on the fine and the coarse grid, and carry the difference found over to the rest (Geological and technical uncertainty of a project).

When upscaling does not work

Sometimes no choice of method and boundary conditions helps: a coarse cell cannot pass on what decides the flow. The signs of such cases, and what is done instead of upscaling, are in the table.

Cases in which a coarse cell does not pass on the flow

Case: No scale separation
Sign: a cell about as large as the correlation length; boundary conditions give different $k^{*}$
What is done: smaller cells, a buffer zone, transmissibility upscaling, global conditions

Case: Thin long barriers and conducting thin layers
Sign: bodies longer than the cell; the $k_v$ of the fine and the coarse grid differ
What is done: explicit barriers, transmissibility multipliers, separate layers

Case: Unstable displacement
Sign: the displacing agent is more mobile than the oil; tongues
What is done: fine grid, a check on a rotated grid

Case: A fractured reservoir
Sign: the fractures conduct, and the matrix stores
What is done: a dual porosity model

Case: Miscible (the fluids mix) or compositional (the composition of the fluids changes) displacement
Sign: the width of the mixing zone depends on the grid
What is done: fine grid, or corrections checked on it

Case: The pattern of wells changes
Sign: conversion to injection, infill drilling
What is done: recalculation for the new pattern, global methods

Case: Cones near wells
Sign: water or gas is drawn up to the wellbore
What is done: local grid refinement, a radial grid near the well

Case	Sign	What is done
No scale separation	a cell about as large as the correlation length; boundary conditions give different $k^{*}$	smaller cells, a buffer zone, transmissibility upscaling, global conditions
Thin long barriers and conducting thin layers	bodies longer than the cell; the $k_v$ of the fine and the coarse grid differ	explicit barriers, transmissibility multipliers, separate layers
Unstable displacement	the displacing agent is more mobile than the oil; tongues	fine grid, a check on a rotated grid
A fractured reservoir	the fractures conduct, and the matrix stores	a dual porosity model
Miscible (the fluids mix) or compositional (the composition of the fluids changes) displacement	the width of the mixing zone depends on the grid	fine grid, or corrections checked on it
The pattern of wells changes	conversion to injection, infill drilling	recalculation for the new pattern, global methods
Cones near wells	water or gas is drawn up to the wellbore	local grid refinement, a radial grid near the well

Unstable displacement and the grid orientation effect are in the lessons Numerical dispersion and Displacement of oil by gas. Fractured reservoirs are in Fractures and naturally fractured reservoirs.

The general way out is to spend cells where they are needed: local grid refinement near wells and at fronts; non-uniform grids along the flow paths; sector models with a fine grid for particular questions; and calculation along streamlines (the paths of the fluid) for a fast run through geological options. Multiscale methods are also being developed. The pressure is solved on a coarse grid, but with functions built from the fine one, and the flows are reconstructed on the fine grid. This is no longer upscaling of properties, but a way to solve the fine problem faster.

Terms and symbols

Scale separation: The condition $\lambda \ll \Delta \ll L$ . The heterogeneity is much finer than the cell, and the flow changes over distances much larger than the cell.
Representative elementary volume (REV): The smallest volume in which an averaged property stops depending on the size of the volume.
Boundary-condition difference: The difference between the values of $k^{*}$ that different boundary conditions give for one cell. It is a sign of cells as large as the correlation length.
Stepwise upscaling (geopseudo functions): The transfer of properties from laminae to a cell of the geological model and onward, by geologically meaningful scales.
Two-point and multipoint flux approximation (TPFA, MPFA): The flow through a face found from the pressures of two cells, or of several. The second kind passes on the full tensor.
K-orthogonal grid: A grid where the vector $\mathbf{K}\mathbf{n}$ of each face points along the line between the centres of the neighbouring cells. In plain words, the flow through each face then depends only on the pressures of the two cells beside it. Only on such a grid is the two-point scheme exact.
Directional permeability $k(\varphi)$: The component of the velocity along a direction at a unit gradient along it: $\mathbf{n}^{\mathsf{T}}\mathbf{K}\mathbf{n}$ .
Well index from a local calculation $\mathrm{WI}^{*}$: The index at which the coarse model gives the rate of the fine one at the same pressures.
Local, global and local–global upscaling: The boundary conditions of a cell calculation come from an assumption (local), from a calculation of the whole model with wells (global), or from an alternation of the two (local–global).
F–Φ curve, first arrival of a tracer: The cumulative share of the flow against the share of the pore volume, from fast paths to slow ones. Its first point gives the arrival of the tracer.

Common mistakes and misconceptions

Wrong

“Cells as large as the correlation length are a reasonable compromise: they are larger than the small ones and smaller than the large ones.”

In fact

It is exactly for them that the equivalent permeability depends most strongly on the boundary conditions. In the calculation of the lesson, at $\sigma_{\ln k}$ = 1.5, for blocks of 1.3 correlation lengths the median difference is 23%. Such cells are checked with two kinds of conditions and with a calculation on a fine grid.

Wrong

“The thinner the layers of the coarse grid, the more accurate the vertical permeability.”

In fact

Each inner face works as a conducting film in the calculation of the $k_v$ of a cell. In the example with lenses 100 m long, the grid of twenty layers of 0.5 m overestimated $k_v$ by 5.0 times, and the one-layer grid was wrong in the median by 2%.

Wrong

“The off-diagonal terms of the tensor are small, so they can be dropped.”

In fact

They are small only if the principal axes are close to the grid axes. In the block of the lesson $k_{xy}$ = 20.0 mD with a diagonal of about 63 mD, but without it the grid loses all the anisotropy: about 63 mD in all directions instead of 83.3 and 43.3.

Wrong

“The well in the model gives less than in the field, so we will choose a rate multiplier in history matching.”

In fact

The correction goes into the well index, by a local calculation on a fine grid. In the example of the lesson the index of a well in a channel is 5.3 times larger than the Peaceman index with the permeability of the cell. A multiplier chosen during history matching hides the cause, and it carries the error into the forecast.

Wrong

“The single-phase rate of the coarse model agreed with the fine one, so the upscaling is right.”

In fact

The rate is a sum. Displacement depends on how the flow is distributed over fast and slow paths. The one-layer model of the eight layers of the lesson passes the same rate, but its tracer arrives after one pore volume, and in the fine model after 0.21.

In depth: advanced detail

Why boundary conditions matter most for cells as large as the correlation length

The boundary conditions influence the flow in a strip near the faces of the block, about one correlation length wide. Farther from a face, the flow is already decided by the rock itself. If the block is much smaller than the correlation length, the rock in it is almost uniform, and the conditions hardly matter. If the block is much larger, the strip takes a small share of its area, about $\lambda/\Delta$ , and the difference falls roughly in inverse proportion to the size of the block. Between these limits the difference passes through a maximum, for blocks of the order of the correlation length (in the calculation of the lesson, the largest of the calculated values is for blocks of 1.3 correlation lengths).

The scatter of $k^{*}$ between blocks also falls with the size of the block. In a plane field, the variance of the average over a block falls roughly in inverse proportion to the number of independent patches in it, that is, as $(\lambda/\Delta)^2$ , and the standard deviation falls as $\lambda/\Delta$ . In the figure, for blocks larger than the correlation length, the share of the scatter that is kept falls quickly: 0.40, 0.23, 0.09.

How the calculation was done. A field of 64 × 64 cells of 25 m is built by the generator of the lesson Reservoir heterogeneity and permeability averaging: a standard normal field with a correlation that falls off exponentially with distance and a decay scale of 2 cells (the correlation length is three scales, 150 m), with $\ln k = \ln 50 + \sigma z$ . The same fields are used for the three values of $\sigma$ . The field is cut into square blocks of 1, 2, 4, 8, 16 and 32 cells. In each block, $k^{*}_x$ is calculated from the two-point cell equations (flow calculation in the block), with sealed side faces and with linear pressure on all the faces. The statistics are over all the blocks of ten realisations. For elongated bodies and three-dimensional flow the numbers will be different. We get them by the same calculation on fields of our own reservoir.

Two-point and multipoint flow through a face

By Darcy's law with a tensor, the flow through the east face of a cell, with the normal along $x$ , is $u_x = -(k_{xx}\,\partial p/\partial x + k_{xy}\,\partial p/\partial y)/\mu$ . The two-point scheme replaces the first derivative by the difference of the pressures of the neighbouring cells, divided by the distance between their centres. It has no term for the second derivative: the pressures of the two cells on either side of the face say nothing about the gradient along the face. If $k_{xy} \neq 0$ , a term is lost that does not get smaller when the cells are made smaller, because the gradient along the face does not vanish. The same happens with a diagonal tensor on a skewed cell: the line between the centres is not perpendicular to the face, and the pressure difference along it contains a component along the face.

Multipoint schemes also restore the gradient along the face, from the pressures of the cells above and below on both sides of the face. In a plane problem the flow through a face then depends on six cells, and the cell equation on the cell itself and its eight neighbours. In three dimensions it is up to twenty-six neighbours. The coefficients are chosen so that the scheme is exact for any linear pressure field and any constant tensor. On a K-orthogonal grid it reduces to the two-point scheme.

Travel time and the F–Φ curve of a model

In the general case, with areal heterogeneity and wells, the F–Φ curve is built from the time of flight: the time in which a particle of fluid travels from the injection well to a given point along a streamline. From the single-phase velocity field we find the time of flight for each cell. We sort the cells by the increasing total travel time from the injection well to the production well. The cumulative share of their pore volume $\Phi$ is plotted against the cumulative share of the rate $F$ that passes through them. For continuous layers without crossflow the travel time of a layer is proportional to $\phi/k$ , and the curve coincides with the Lorenz curve of the layers.

This diagnosis of the flow needs only one pressure calculation and one set of streamlines. So we do it on the fine and the coarse model, and we compare the curves, the Lorenz coefficient and the time of first arrival at each production well. The difference shows which fast paths the coarse grid has lost, even before the displacement calculations.

Connections to other topics

Permeability upscaling for single-phase flow — equivalent permeability, boundary conditions, transmissibility upscaling and a first estimate of the well index.
Two-phase upscaling: pseudo-functions — pseudo functions, their methods and limitations.
Reservoir heterogeneity and permeability averaging — correlation length, anisotropy, the Lorenz curve.
Reservoir architecture and performance — connectivity of bodies, barriers and baffles, types of architecture.
Well model in a grid block — the well index and the equivalent radius of a cell.
Grids in reservoir simulation models — corner-point grids, local grid refinement, the grid orientation effect, convergence.
Petrophysical input for simulation: relative permeability and capillary pressure — relative permeability and capillary pressure curves of laminae for the first step.
Continental environments: aeolian and fluvial — inclined laminae and the permeability tensor of a cross-bedded set.
Pore-scale flow modelling — the representative volume and the properties at the pore scale.
History matching — why upscaling errors must not be hidden in multipliers.

Check yourself

1
A cell of 200 m, a correlation length of 150 m, $\sigma_{\ln k}$ ≈ 2. The flow calculation in the cell gave 40 mD with sealed side faces and 70 mD with linear conditions. What does this mean?
1. A.We should take the average, 55 mD, which is more accurate than either
2. B.We should take 70 mD, because linear conditions are always more accurate
3. C.A calculation error: the right conditions give one value
4. D.The cell is about as large as the correlation length, and its equivalent permeability is not defined uniquely. We need a buffer zone, a different cell size or a check on a fine grid
Show the answer
Correct answer: D. The cell is about as large as the correlation length, and its equivalent permeability is not defined uniquely. We need a buffer zone, a different cell size or a check on a fine grid
The ratio of the cell size to the correlation length is about 1.3, exactly where the difference between the boundary conditions is largest. At $\sigma_{\ln k}$ ≈ 2 the median difference in the calculation of the lesson is about a half. The value depends on whether the bodies continue in the neighbouring cells, and the cell does not know this. The average is not more accurate: the true value may lie closer to either end. Linear conditions give the largest value, but not necessarily the right one.
2
A unit with discontinuous shale lenses was coarsened horizontally, keeping layers of 0.5 m. The vertical permeability of the coarse model came out five times larger than that of the fine one. Why?
1. A.The sealed side faces of the cells underestimate $k_v$
2. B.In the calculation of the $k_v$ of each cell, the pressure on its top and bottom faces is uniform, like a conducting film along which the flow freely goes round the lenses
3. C.Numerical dispersion is to blame
4. D.Thin layers pass the lenses better, so the fine model is wrong
Show the answer
Correct answer: B. In the calculation of the $k_v$ of each cell, the pressure on its top and bottom faces is uniform, like a conducting film along which the flow freely goes round the lenses
In the fine model a flow that has passed a gap in one lens must still reach the gap in the next one, and the path gets longer. In the coarse model the vertical permeability of each cell is calculated with the same pressure over the whole width of its faces, and the cells add up in series. There are no long detours, and the result approaches the rule “share of sandstone × permeability”. Sealed side faces act in the opposite direction and underestimate $k_v$ , but in thin cells their role is small. Numerical dispersion belongs to the transport of saturation, not to single-phase permeability.
3
The principal permeabilities of a cross-bedded set are 400 and 40 mD, and the principal axes are turned by 45° to the grid axes. Find $k_{xx}$ , $k_{yy}$ , $k_{xy}$ and the directional permeability along the laminae that the two-point scheme will see.
Show the answer
At 45° $\cos^2 = \sin^2 = 0.5$ and $\sin\cos = 0.5$ . So $k_{xx} = k_{yy} = 0.5\cdot400 + 0.5\cdot40 = 220$ mD, and $k_{xy} = (400 - 40)\cdot0.5 = 180$ mD. A scheme without $k_{xy}$ sees $220\cdot0.5 + 220\cdot0.5 = 220$ mD in any direction. Along the laminae it is 45% less than 400 mD, and across them 5.5 times more than 40 mD. The tenfold anisotropy disappears completely.
4
A local calculation near a well on a region of 3 × 3 cells gave an index of 410 mD·m, and a ratio of the rates of the fine and the coarse models of 2.6. On a region of 5 × 5 it gave 418 mD·m and 2.2. What should be carried into the full model?
1. A.An index of about 415 mD·m: it hardly depends on the region, while the ratio of the rates includes the resistance of the distant cells
2. B.Nothing: the results for different regions contradict each other
3. C.A rate multiplier of 2.2, because the region is larger
4. D.A rate multiplier of 2.6, found on the smaller and more accurate region
Show the answer
Correct answer: A. An index of about 415 mD·m: it hardly depends on the region, while the ratio of the rates includes the resistance of the distant cells
The index describes the neighbourhood of the well and changes by 2% with the region. The ratio of the rates depends on how much of the total resistance comes from the distant cells, which the coarse model passes on correctly. So it changes with the size of the region, and it is not carried into the full model.
5
A unit of two continuous layers: 2 m, 500 mD, $\phi$ = 0.25 and 6 m, 50 mD, $\phi$ = 0.20. When will the tracer arrive in single-phase injection through the fast layer if the layers are not combined, and if they are combined into one layer of the coarse grid?
Show the answer
The sum $\phi h = 0.5 + 1.2 = 1.7$ m, and the sum $kh = 1000 + 300 = 1300$ mD·m. The fast layer is the first one ( $k/\phi$ = 2000 against 250 mD): $\Delta C = 0.5/1.7 = 0.294$ , $\Delta F = 1000/1300 = 0.769$ , and the arrival is after $0.294/0.769 \approx 0.38$ pore volumes. In the combined layer the tracer arrives after 1 pore volume, which is 2.6 times later. The single-phase rate is the same: the $k$ of the combined layer is the arithmetic mean, 162.5 mD.
6
An injection well takes 200 m³/d into a reservoir 10 m thick. At a Darcy velocity of 0.02 m/d the ratio of gravity to viscous forces of a coarse cell is 0.4. At what distance from the well will it reach one?
Show the answer
$N_{gv}$ is inversely proportional to the velocity. It is one at $u = 0.02\cdot0.4 = 0.008$ m/d. By the formula $u = q/(2\pi rh)$ , $r = 200/(2\pi\cdot10\cdot0.008) \approx 398$ m. Closer to the well the cell is in the viscous regime. Farther away gravity has time to separate the phases, and the pseudo functions of these zones must differ.

Updated: 19 October 2026