Simulation of contaminant transport in an oil-bearing stratum at water flooding

SIMULATION OF CONTAMINANT TRANSPORT IN AN

OIL-BEARING STRATUM AT WATER FLOODING

1

Boris N. Chetverushkin *, Natalia G. Churbanova * _{, Anton A. Soukhinov}† _{and Marina} A. Trapeznikova *

*Institute for Mathematical Modeling RAS 4-A Miusskaya Square, Moscow 125047, Russia

e-mail: nata@imamod.ru

† _{Moscow Institute of Physics and Technology}

9, Institutskii per., Dolgoprudny, Moscow Region, 141700,Russia

e-mail: Soukhinov@mail.ru

Key words: Porous Media Flow, Multiphase Contaminant Transport, Adaptive Cartesian Mesh

Abstract. The research deals with simulation of multiphase immiscible incompressible fluid flows in porous media when one of liquids contains a contaminant. The model neglecting capillary and gravity forces (the Buckley-Leverett model) is used. The problem of passive contaminant transport in an oil-bearing stratum at water flooding is investigated. The numerical algorithm simulating flows in heterogeneous porous media is developed. The procedure of adaptive Cartesian mesh generation is constructed.

1 INTRODUCTION

The work is aimed at development of efficient algorithms for simulation of multiphase fluid flows in porous media when one of fluids contains a contaminant. Simulation of such flows is of great practical importance because it is necessary to predict them, for example, while developing oil recovery technologies or investigating ecology problems. As a typical example the problem of passive contaminant transport in an oil-bearing stratum at water flooding is treated. It is supposed that the oil field is covered by a network of water injection and oil production wells and the contaminant arrives with water through a number of wells (for example, salty water is injected at some sources).

Numerical implementations are based on finite difference approximations on rectangular grids of the MAC type. Symmetric and periodic boundary conditions are implemented. Implicit algorithms developed by the authors earlier3, 4 are mainly used. A natural extension of them for heterogeneous media is presented here. The above problem was solved as for homogeneous as for heterogeneous fields with different well configurations. The accuracy was estimated by comparing results obtained on computational grids of different coarseness and by checking the oil balance.

A widely used tool for the accuracy increase is employment of the adaptive mesh refinement5-7. One of objects of the present work was development of original algorithms for the embedded mesh generation with adaptation to the problem solution by splitting and merging cells.

2 PASSIVE CONTAMINANT TRANSPORT IN AN OIL-BEARING STRATUM At the present time exploitation of most oil and gas deposits are at the closing stage and in practice secondary methods, in particular the water flooding, are mainly used for oil and gas recovery. Numerical simulation of these technological methods is usually aimed at developing an optimal strategy of oil reservoir exploitation.

Let us consider a joint flow of oil and water in an oil field when oil is extracted by means of non-piston water displacement1, 2. In this problem formulation an oil-bearing stratum is assumed to be thin so that the 2D (plane) problem is considered. The oil reservoir is covered by a regular network of vertical wells, namely, water injection wells (sources) and oil production wells (sinks). Their disposition schemes may be quite different. Reference lengths in the problem of oil recovery are hundreds and thousands of meters and reference times are months and years. In these conditions the ratio of the capillary pressure to the full hydrodynamic pressure drop is small. It allows to omit capillary forces so that the flow of two phases w (water) and o (oil) in a porous medium is governed by the classical Buckley-Leverett model based on a number of assumptions: fluids are immiscible and incompressible, the medium is undeformable, phase flows comply with the Darcy law, capillary and gravity forces are neglected.

For computations the complete mathematical statement is reduced to the next system of equations:

( ) ( )

(

)

_{( )}

( )

at sources

div

( )

( )

( )

( )

w w w w w w o w w o k s F s k s k s µ µ µ = + (3)

( )

w

( )

w o

( )

w w w o k s k s K s k µ µ   = − _ + _   (4)

The water saturation sw and the pressure P are sought in a subdomain of symmetry cut

from an unbounded uniform stratum. The next notations are used here: m is the porosity, q describes debits of wells, F(sw) is the Buckley-Leverett function, K(sw) is a nonlinear

coefficient including the absolute permeability k , relative phase permeabilities kw(sw) , ko(sw)

and dynamic viscosities µw , µo . Relative permeabilities are defined by the next relations

( s = 0.1 - the bound saturation, s = 0.8 - the critical saturation, n1 = 2, n2 =2, n3 = 0.5 in

our tests):

( )

1

( )

2 3 1 1 1, 0 0, 0 , , 0, 1 0.8 , w w n n w w o w w w w w n w _w w s s s s s s s s k s s s s k s s s s s s s s s s _{s s} s s s s s   ≤ < ≤ <  _  _  −   −    =_{ } ≤ ≤ =_{  } ≤ ≤ − −        _{< ≤}  _{ − }  _{ × < ≤}    _ − _  (5)

All numerical implementations are based on finite difference approximations on rectangular grids of the MAC type.

Solution of transport equation (1) faces significant difficulties due to existence of a discontinuity in the saturation function. Different explicit as well as implicit numerical methods for solving this equation were studied3. In the current paper equation (1) is solved by the iterative secant method what demonstrates advantages over explicit difference schemes.

Equation (2) in our implementation can be solved by the original iterative

(

α β−

)

-algorithm or by the local relaxation (ad hoc SOR) method with red-black data ordering4.

(

α β−

)

-algorithm is an extension of the Tri-Diagonal Matrix Algorithm to the two dimensional case. To reduce the number of iterations, the successive overrelaxation is applied to this algorithm and so-called

(

α β−

)

with relaxation is obtained.

flow in a porous medium is governed by the same equations (1)-(2) and the Darcy velocity of water Ww is obtained via the saturation and the pressure (the generalized Darcy Law):

( )

= - w w w w k s k P µ W grad (6)

Adsorption and convective diffusion have mainly an influence on contaminant spreading in porous media. Then the equation for the contaminant concentration can be written like this:

(

w

)

_div

₍

_{( )}

₎

w w c s c a m c Q t ∂ + + + = ∂ W S W (7)

Here c is the contaminant solution concentration, a is the concentration of the contaminant being adsorbed in a unit volume of a porous medium (a is a function of c) so that a unit volume of a porous medium filled by fluid with a contaminant contains (mc + a) moles of the contaminant. Qc describes a source of the contaminant. S is a “diffusion flow” caused by the

convective diffusion:

{ }

, , , i i j j c S D i j x y x ∂ = − ∈ ∂ (8)

Dij is an effective tensor of the convective diffusion. There are phenomenological

formulas to obtain it via the Darcy velocity and some characteristics of the medium. In our computations we use the formula after V.N.Nikolaevskij8:

(

1 2

)

2

i j i j i j w

D =_ λ λ δ− +λ n n W_ (9)

where W_w = W_w , ni is the unit vector in the direction of the water Darcy velocity, λ1 and λ2

are some positive coefficients, their order of magnitude agrees with the character size of the medium microheterogeneity. One can see that even in an isotropic medium there is a direction singled out by the Darcy velocity vector.

Adsorption is not taken into account yet. For discretization of equation (7) the finite difference scheme with central differences is used. As equation (7) contains mixed derivatives in term div(S) (with the account of (8), (9)) the nine-point stencil is used for approximation.

The computational domain is closed to a tore in the both directions in order to provide circulating boundary conditions. The two-level scheme with the weight 1/2 is proposed to provide the second order of approximation on time as well as on space:

(

) (

)

(

)

1 1 , , , , 1 _{, ,} 1_, 1_, 1 2 n n n n i j i j i j i j n n n n n n s c s c m F c W f F c W f t + + + + + − = + ∆ (10)

Since the history of the flow is known values depending on Ww is taken not only from the

Equation (10) is solved by the local relaxation (ad hoc SOR) method with special data ordering. As the nine-point stencil is used the red-black ordering is not suitable in this case and multicolor data ordering is employed.

As test problems different variants of well disposition were considered. The figures below illustrate the so-called honeycomb well configuration. At the initial moment the water saturation equals the bound saturation s, the pressure is constant all over the domain, the concentration equals zero. A passive contaminant arrives into the homogeneous stratum with water through the central injection well. Figure 1 presents pressure, water saturation and concentration fields for different time moments.

Pressure Water saturation Concentration

(a) –the 10th _day

Pressure Water saturation Concentration

Pressure Water saturation Concentration

(c) - 100th _day

Pressure Water saturation Concentration

(d) - 200th _day

Figure 1: Pressure, water saturation and concentration fields at different time moments for the honeycomb well disposition

Figures 1 (a), (b) reflect beginning stages of the process. Contaminant spreading is some behind of water saturation spreading. Figure 1 (c) demonstrates that the concentration reaches production wells. And Figure 1 (d) shows that the concentration can not spread further under such initial conditions and at such specified debits of the contaminant on injection wells and debits of the mixed fluid on production wells.

3 HETEROGENEOUS POROUS MEDIA

considered as a function of coordinates but not as a tensor that is the anisotropy is not taken into account yet.

Several test problems were solved in heterogeneous media. In any case a contaminant arrives into the stratum with water. At first the computational domain containing one water injection and one oil production wells was considered with the circulating boundary conditions. In computations these conditions are provided by closing the domain to a tore in the both directions. Then the classical five-spot pattern with the boundary conditions of symmetry was considered. The model configuration consists of a square with a water injection well in the center and oil production wells in corners. The symmetry allows for the reduction of the problem to the right upper quadrant of the domain. Figure 2 shows the permeability field, the water saturation and the concentration plots at some time moment. The maximum saturation (the red color) corresponds to the injection well. One can see that the developing fronts have a complex shape and the so-called “fingering” takes place.

Figure 2: The permeability, the water saturation and the contaminant concentration for the five-spot problem.

In order to estimate the accuracy of the developed algorithms results obtained on computational grids of different coarseness are compared. Figure 3 presents the water saturation fields over grids of 100 × 100 and 200 × 200 nodes. A good agreement is achieved.

Predictions for the seven-spot (“honeycomb”) well configuration was also performed. A contaminant arrives only through the central water injection well. The corresponding water saturation and the contaminant concentration in the given heterogeneous stratum at some time moment are depicted in Figure 4 .

Figure 4: The permeability, the water saturation and the contaminant concentration for the seven-spot problem.

4 RECTANGULAR ADAPTIVE MESH FOR SOLVING CONVECTION AND FILTRATION PROBLEMS

4.1 Two-dimensional adaptive mesh

Let's consider a 2D adaptive mesh and algorithms of its functioning. Cells are rectangular. We shall allow each cell to be divided into four cells of the identical size. Four adjacent cells of the identical size can be combined in one. For simplification of algorithm we allow merging only those cells which once made one cell. We start splitting from one square cell

0

C , which is referred to as a root of a mesh tree. Each cell C (_i

0

≤

i

≤

n

−

1

) stores real number N , describing mean value of calculated function within the cell. _i

The mesh can be stored in the form of quaternary tree in which nodes are cells which have been split, and leaves are cells on which calculations on a current time layer are made. An example of such a tree is shown in the Figure 5:

When the cell C of the mesh is split into 4 child cells _i C_LUi,C_RUi,C_LBi,C_RBi, they get values that equality

(

N_LU N_RU N_LB N_RB

)

N_i

i i i

i + + + 4= takes place. When cells are merged,

they are destroyed, and their average value is stored to their parent. It is important in order that repeated splitting/merging of cells did not lead to error accumulation.

However it is not enough to be able to merge and split cells. It is still necessary to solve following problems:

• It is necessary to determine what cells are around of the current cell.

• It is required to determine values, which need to be attributed to cells at splitting to most precisely solve the problem.

• It is needed to approximate differential equation on such a mesh.

• It is required to define, which cells to agglomerate and which to break to obtain precise solution.

The first problem is resolved by means of the algorithm of neighbors determination. Two following problems are solved via interpolation whereas the last problem is solved through the refining algorithm.

4.2 Neighbors determination algorithm

To interpolate and refine the mesh, we will need to find neighbors for each cell.

Let the cell of the mesh C stores indexes _i

LU

_i

,

RU

_i

,

LB

_i

,

RB

_i of the child cells and the index

P

_i to the parent cell. If some of the listed cells are absent, corresponding indexes are negative. Then for finding of neighbor cells it is suggested to use the following recursive formula (example of finding of the left neighbor for the case of periodic boundary conditions):

( )

( ) ( )

( )

_{( )} ( ) ( )

( )

_{( )}



























=

≥

<

=

≥

≥

=

≥

=

≥

<

=

≥

≥

=

≥

<

=

.

,

0

if

,

;

0

,

,

0

if

,

;

0

,

,

0

if

,

;

,

0

if

,

;

0

,

,

0

if

,

;

0

,

,

0

if

,

;

neighbor)

a

as

itself

has

tree

the

of

(root

0

if

,

i

RB

P

LB

RB

i

LB

P

P

LEFT

RB

i

LB

P

RB

i

RU

P

LU

RU

i

LU

P

P

LEFT

RU

i

LU

P

RU

P

i

i

LEFT

i i i i i i i i i i i i i i P i P P LEFT P i i P LEFT P i P LEFT P i P P LEFT P i i P LEFT P i P LEFT i (11)

Figure 6: An example of work of neighbor finding algorithm

4.3 Interpolation

Interpolation is applied to leaves of the mesh tree. Using values of the mesh cells we can construct a continuous function using interpolation. This continuous function should satisfy the following conditions: its mean values on leaves C should be equal to the values _i N , _i which stored in these cells.

First of all, we will restrict sizes of cells of the mesh: sizes of cells having common points should not differ more than twice.

For interpolation we split each cell C into four square bilinear patches which are based _i on nine points _pi _pi

8

0,..., (figure 5). Values Pji in points pij are calculated on the basis of values i

j

n

N of the neighbor cells i j

n

C numbered i j

n . We will omit index i of the current cell in this section. Instead of i

j

n

N we will write N . The cell C_j can have from 6 to 12 neighbor cells (Figure 7): p1 p2 p3 p0 p4 p5 p6 p7 p8 n0 n1 n2 n3 n4 n5 n6 n7 n8 n9 n10 n11 n12

If the cell has less than 12 neighbors some of numbers n will be equal. The neighbour _j cells n₂,n₅,n₈,n₁₁ we shall name angular neighbors of the cell n , other neighbor cells - ₀ boundary neighbors.

Let's result reasons which will allow us to calculate values P_j, j=0,...,8.

First of all, interpolation is applied to calculation of new cells values at cells splitting. Their value is equal to average value of interpolation function on the cell, therefore values in four new cells are calculated using the following formulas:

4 0 5 7 1 P P P P N_LU = + + + , 4 6 0 2 7 P P P P N_RU = + + + (12) 4 8 3 0 5 P P P P N_LB = + + + , 4 4 8 6 0 P P P P N_RB = + + + . (13)

Then to not conflict with average value on the cell, values P should satisfy the following _j restriction:

(

)

16 2 4 4 8 7 6 5 4 3 2 1 0 0 P P P P P P P P P N N N N N ₌ LU + RU + LB+ RB ₌ + + + + + + + + (14)

Here N is the value in the current cell. Therefore, if values ₀ P₁,...,P₈ will be known, we can calculate value P : ₀

(

)

4 2 4 1 2 3 4 5 6 7 8 0 0 P P P P P P P P N P = − + + + + + + + (15)

Further, cells n , ₀ n₁, n₂ and n , being around of point₃ p₁, possess common information only about values N₀,N₁,N₂,N₃. Let's not consider the sizes of cells. Then that interpolation function was continuous, we can take average value of adjoining cells. The same reasoning concerns other three angular cells. Therefore following formulas are used:

. 4 , 4 , 4 , 4 12 11 10 0 4 9 8 7 0 3 6 5 4 0 2 3 2 1 0 1 N N N N P N N N N P N N N N P N N N N P + + + = + + + = + + + = + + + = (16)

If the cell has two neighbors along a border (for example n₁ ≠n₇), then values on border of the cell should coincided with values in corners of the neighbor cells. We should calculate

8 5,...,P

. 4 2 ˆ , 4 2 ˆ , 4 2 ˆ , 4 2 ˆ 0 10 9 8 0 4 3 7 0 12 6 6 0 7 1 5 N N N P N N N P N N N P N N N P + + = + + = + + = + + = (17)

If boundary neighbor is larger than current cell (for example, n₁ =n₇ =n₈ or n₂ =n₁ =n₇), for ensuring of interpolation function continuity we should use following values:

. 2 , 2 , 2 , 2 4 3 8 2 1 7 4 2 6 3 1 5 P P P P P P P P P P P P = + = + = + = + (18)

It is necessary to consider a case when boundary neighbor has the same size as current cell. Applying one-dimensional analogue of the described interpolation to boundary edges, we will get: 8 6 3 3 ~ , 8 6 3 3 ~ 8 6 3 3 ~ , 8 6 3 3 ~ 10 9 12 11 8 7 8 0 4 3 6 5 2 1 7 12 6 11 10 5 4 6 0 7 1 9 8 3 2 5 N N N N N N N P N N N N N N N P N N N N N N N P N N N N N N N P + + + − − − − = + + + − − − − = + + + − − − − = + + + − − − − = (19)

Then the right values for P₅,...,P₈ can be chosen in the following way (on example of P ₅ value):      = = = = ≠ = . otherwise , ~,if or ; ; if , ˆ 5 7 1 2 8 7 1 5 7 1 5 5 P n n n n n n P n n P P (20)

Interpolation is needed not only to finding values of new cells at splitting, but also to do approximation of the differential equation and in mesh refining algorithm. We will describe briefly a way of approximation of the differential equation. The equation can be approximated by a method of finite differences: corresponding continuous derivatives are replaced by their difference analogues, for example:

( )

( )

( )

₍

₎

( )

, 1 _. , 2 2 , , , , , 3 4 1 2 2 2 6 0 5 2 2 5 6 0       − ₋ − → ∂ ∂ ∂ + − → ∂ ∂ − → ∂ ∂ → x x y x x h P P h P P h x y y x U h P P P x y x U h P P x y x U N y x U (21)

Besides that interpolation can be used for transition from one mesh to another. It is possible to store various physical values in various adaptive meshes. Then, if for calculation of value in a cell C from mesh MESH₁ value from other mesh MESH₂ is needed, it is possible to numerically integrate interpolation function of mesh MESH₂ over the cell C. Cells of mesh MESH₂, which have got in C entirely, will enter into integral proportionally to their area. And the cells crossed with C partially, should be interpolated, and integrals from bilinear patches are taken proportionally to their areas.

4.4 Mesh refinement algorithm

To employ successfully an adaptive mesh, we should use computationally simple algorithm of the mesh adaptation (“refining”) to the computed solution. To describe algorithm of mesh refinement, the concept of data variation of a cell is required.

Let's name data variation in the cell C following value: _i

(

1

)

₂ , 0 2 0 2 2 0 0 x y yi xi i i i h h h h m M m M + + + − − − = ∆ ε ε (22) Here

(

)

    = = ( ),if cellwhole. max ; dividual cell if , , , , max 8 ,..., 0 i j j RB LB RU LU i _P M M M M M i i i i

(

)

    =

=min ( ),if cellwhole.

; dividual cell if , , , , min 8 ,..., 0 i j j RB LB RU LU i _P m m m m m i i i i i i i i RU LB RB

LU , , , - numbers of child cells of cell C . _i

yi xi h

h , - sizes of cell C . _i

i j

P - interpolation points of cell C . _i

<

ε - the parameter describing influence of the cell size on necessity of its splitting. The cell which number is zero is the root of the mesh tree.

Let's admit that mesh refinement algorithm decided to break some cell. But this cell at the moment of splitting can adjoin to larger cells; that does splitting impossible because adjoining cells cannot differ in sizes more than twice. In this connection it is necessary to split in the beginning large neighbors of the cell which, in turn, can lead splitting of another cells. We shall speak therefore about the set of split cells which generated by splitting of a cell C . _i

be realized by merging of unnecessary cells: the cell with small variation having child cells can be united by destroying the child cells. We cannot unite each cell: in its neighborhood can be small cells which need to be merged before uniting of the initial cell and so on. However there is no necessity to make a chain of merges: it is enough to try to unite last element of this chain. If it is impossible because of big variation, then the chain especially cannot be merged. We will say that cell

C

_i can be united if it doesn’t has neighbors smaller than

C

_i.

So, we offer algorithm which reduces the maximal variation of the mesh by splitting chains of the cells containing the cell with maximal variation. For maintenance of necessary quantity of cells, cells having the minimal variation are merged. The algorithm stops as soon as the minimal variation of the cell to unite becomes more than variation of the cell which is a subject to splitting. In the last case no action can reduce the maximal variation of data in the mesh any more.

n

is the number of mesh cells,

n

is the needed number of mesh cells. Here is the algorithm:

1. While

n

>

n

and we can unite cells we are uniting cells having minimal variation; if we can’t unite cells and still

n

>

n

the algorithm stops. When

n

≤

n

we can perform refining (the next steps of the algorithm).

2. We find cell

i

_split with maximum variation and determine the number

n

_split of new cells in the split chain.

3. If

n

+

n

_split

≤

n

we split cell

i

_split and its chain, and go to step 2.

4. If

n

+

n

_split

>

n

we find cell

i

_merge with minimal variation which can be united. If we have found such a cell and

split

merge i

i

≤

∆

∆

, we merge this cell. Otherwise algorithm stops.

5. The merging of cells can not affect the cell

i

_split, so we do not need to find

i

_split again. Go to step 3.

4.5 Numerical results on adaptive meshes

Let's consider an example of calculation of a convection problem on homogeneous and adaptive meshes. The calculation domain is the unit square. The input condition is a spot in the form of circle in diameter 0.25 which center has coordinates

(

0.25,0.25

)

. The velocity is vector

( )

1,1 in all points. At time 0.7 spot should reach coordinate

(

0.75,0.75

)

.

X Y 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X Y 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 8: The entry condition for homogeneous and adaptive meshes

X Y 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 X Y 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Figure 9: The result of calculation on homogeneous and adaptive meshes

5 CONCLUSIONS AND FUTURE WORK

– The problem of passive contaminant transport in an oil-bearing stratum at water flooding in a homogeneous and heterogeneous media has been studied numerically. – A new algorithm of 2D adaptive mesh refinement was constructed, validated on a

– The algorithm of adaptive mesh refinement will be implemented in the form of package in order to use it widely in applications.

REFERENCES

[1] R.E. Ewing (Ed.), The mathematics of reservoir simulations, SIAM, Philadelphia, (1983). [2] V.M. Entov and A.F. Zazovskij, Hydrodynamics of processes for oil production

acceleration, Moscow, Nedra, (1989). - In Russian.

[3] M.A. Trapeznikova and N.G. Churbanova, “Simulation of multi-phase fluid filtration on parallel computers with distributed memory”, In: CFD’98, Proc. of the 4th European CFD Conf. (Eds K.D. Papailiou et al.), Wiley, Chichester, Vol. 1, Pt. 2, pp. 929-934, (1998).

[4] M.A. Trapeznikova, N.G. Churbanova and B.N. Chetverushkin, “Parallel elliptic solvers and their application to oil recovery simulation”, In: HPC‘2000, Grand Challenges in Computer Simulation, Proc. of Simulation Multi Conf., SCS, San Diego, CA, pp. 213-218, (2000).

[5] M.J. Berger and M.J. Aftosmus, “Aspects (and Aspect Ratios) of Cartesian Mesh Methods“, In: Proc. of the 16th Int. Conf. on Numerical Methods in Fluid Dynamics, (6-10 July, 1998, Arcachon, France), Springer-Verlag, Heidelberg, Germany.

[6] J.A. Trangenstein, “Multi-scale iterative techniques and adaptive mesh refinement for flow in porous media”, Advances in Water Resources, 25, 1175-1213 (2002).

[7] L.H. Howell and J.B. Bell, “An adaptive mesh projection method for viscous incompressible flow”, SIAM J. Sci. Comput. 18, N4, 996-1013 (1997).