You must log in to edit PetroWiki. Help with editing

Content of PetroWiki is intended for personal use only and to supplement, not replace, engineering judgment. SPE disclaims any and all liability for your use of such content. More information

# Diagonalizing the permeability tensor

The form of Darcy’s law that is most widely used in formulating fluid flow equations in reservoir simulators assumes that the coordinate system is aligned with the principal axes of the permeability tensor. The resulting diagonalized permeability greatly simplifies the fluid flow equations. The simplified equations are easier to code and can be solved with less computation time than fluid flow equations that include the full permeability tensor.

Research in naturally-fractured-reservoir modeling,^{[1]} geomechanics,^{[2]} and upscaling^{[3]} has demonstrated that the full permeability tensor is needed to correctly solve fluid flow problems in a variety of realistic settings. The mathematical procedure for diagonalizing the permeability tensor is presented here as an illustration of the mathematics discussed in matrices and linear algebra. The relationship between the diagonalized-permeability-tensor assumption and grid orientation is discussed later on this page. An understanding of the relationship between grid orientation and the permeability tensor can help us decide how to orient a fluid flow grid to most accurately represent the permeability distribution in a reservoir. The directional dependence of permeability and the permeability tensor are first introduced. The procedure for diagonalizing the permeability tensor is then presented.

## Contents

- 1 Darcy’s law and the permeability tensor
- 2 Similarity transformations
- 3 Matrix diagonalization procedure
- 4 Diagonalizing a symmetric 2 x 2 matrix
- 5 Eigenvectors
- 6 Coordinate transformation
- 7 Rotational transformation of a 2 x 2 permeability tensor
- 8 Nomenclature
- 9 References
- 10 Noteworthy papers in OnePetro
- 11 Noteworthy books
- 12 External links
- 13 See also

## Darcy’s law and the permeability tensor

In one dimension, Darcy’s law says that flow rate is proportional to pressure gradient. This can be expressed in oilfield units for single-phase flow as

where is flow rate (B/D), *x* is length (ft), A is cross-sectional area (ft^{2}), *μ* is fluid viscosity (cp), *k* is permeability (md), and Φ is the phase potential (psia).

In **Eq. 2**, Δ*z* is depth from a datum (ft), *P* is fluid pressure (psia), and *γ* is the pressure gradient associated with the gravity term (psia/ft). The form of Darcy’s law with full permeability tensor in Cartesian coordinate system {*x*, *y*, *z*} is

where we have treated the cross-sectional area, *A*, as a constant with respect to direction. **Eq. 3** can be rewritten as either a dyadic equation,

by treating permeability as a dyadic, or as a matrix equation,

by treating permeability as a matrix. We are interested here in the matrix representation.

The diagonal permeability elements {*k*_{xx}, *k*_{yy}, *k*_{zz}} represent the dependence of flow rate in one direction on pressure differences in the same direction. The off-diagonal permeability elements {*k*_{xy}, *k*_{xz}, *k*_{yx}, *k*_{yz}, *k*_{zx}, *k*_{zy}} account for the dependence of flow rate in one direction on pressure differences in orthogonal directions. Expanding **Eq. 3** gives the corresponding set of three equations demonstrating this dependence.

## Similarity transformations

**Eq. 5** relates flow rate and the pressure gradient term, . We can use a similarity transformation to diagonalize the matrix while preserving the form of the relationship between and . Let us first show that a similarity transformation preserves the form of **Eq. 5**.

We begin by multiplying **Eq. 5** on the left by to find

where is a nonsingular, *n* × *n* square matrix. Because is nonsingular, it is invertible; that is, it satisfies the equality

where is the *n* × *n* identity matrix. Substituting **Eq. 8** into **Eq. 7** gives

Defining the transformed matrices

and using the similarity transformation

in **Eq. 9** yields

**Eq. 12** is the same form as **Eq. 5**.

## Matrix diagonalization procedure

It is mathematically possible to find a coordinate system {*x*′, *y*′, *z*′} in which the permeability tensor has the diagonal form . We diagonalize the matrix by finding and applying a similarity transformation matrix . The procedure for finding a matrix that diagonalizes an *n* × *n* matrix is as follows:^{[4]}

- Find the eigenvalues {
*λ*_{i}:*i*= 1, ...,*n*} of from the eigenvalue equation det = 0. - Find
*n*linearly independent eigenvectors . - Form the similarity transformation matrix with the eigenvectors as columns.
- Calculate the diagonalized matrix ´. The diagonal entries of ´ are the eigenvalues corresponding to the eigenvectors .

The coordinate axes {*x*′, *y*′, *z*′} are the principal axes of the diagonalized tensor, and the diagonal form of the permeability tensor is obtained by a principal-axis transformation. The flow equations along the principal axes are

The form of the permeability tensor depends on the properties of the porous medium. If the medium is anisotropic, at least two elements of the diagonalized permeability tensor are not equal. If permeability does not depend on direction, then permeability is isotropic, and the elements of the diagonalized permeability tensor are equal, that is,

If the magnitude of the elements of the permeability tensor varies from one point in the medium to another, the permeability tensor is heterogeneous; otherwise, permeability is homogeneous. The principal axes of the permeability tensor may also vary from point to point in the medium if permeability is heterogeneous.

## Diagonalizing a symmetric 2 x 2 matrix

The ideas previously presented are implemented by applying the matrix diagonalization algorithm to the 2 × 2 symmetric matrix

as viewed in the two-dimensional (2D) Cartesian coordinate system ={ x_{1},x_{2}} shown in **Fig. 1**. For this example, we require that the elements of satisfy the relations

The relation *k*_{12} = *k*_{21} for off-diagonal elements is necessary to assure that the matrix is symmetric. The requirement that is symmetric is important when we consider a coordinate transformation. To find the diagonal matrix ´ corresponding to , we must first solve the eigenvalue problem

The two characteristic roots or eigenvalues *λ*_{+} and *λ*_{–} of **Eq. 17** are the diagonal elements of the diagonalized 2 × 2 matrix ´, thus

where *λ*_{+} and *λ*_{–} are calculated as the solutions to the eigenvalue problem.

The eigenvalue problem is det = 0. Using **Eq. 15** gives

or

We expand the characteristic equation to get

The two eigenvalues are found from the quadratic equation to be

The sum of the eigenvalues satisfies the relation

## Eigenvectors

The matrix is composed of orthonormal eigenvectors found from

with the identity matrix . Expanding **Eq. 25** gives

for the eigenvalue *λ*_{+}, and

for the eigenvalue *λ*_{–}. Rearranging **Eq. 26** gives

**Eq. 28** and the normalization condition,

provide the two equations that are necessary for determining the components of ^{+}; thus,

and

Similar calculations for ^{-} yield the results

and

where the relation

from **Eq. 33**, has been used.

To show that ^{+} and ^{-} are orthogonal, we must show that

Substituting in the expressions for ^{+} and ^{-} gives

as expected.

## Coordinate transformation

We now use the orthonormal eigenvectors ^{+} and ^{-} to construct the transformation matrix . According to the algorithm for diagonalizing a square matrix presented previously, we form as

or

A coordinate vector in the transformed coordinate system is given by . Rewriting the matrix equation for coordinate transformations in algebraic form gives

or

An angle (*θ*) can be associated with the linear transformation by writing the 2D coordinate transformation as

The coordinate systems and are related by the counterclockwise rotation shown in **Fig. 1**.

Equating elements of the transformation matrix in **Eqs. 40** and **41** gives

and

The equalities *a*_{1}^{+} = *a*_{2}^{–} and *a*_{2}^{+} = –*a*_{1}^{–} are in agreement with **Eqs. 30** through **33**. The angle *θ* may be found from either

or

Note that the equality,

demonstrates that the eigenvectors are orthonormal.

## Rotational transformation of a 2 x 2 permeability tensor

We want to calculate changes to the permeability tensor when we transform from a coordinate system where only the diagonal elements of a square matrix ´ are nonzero to a coordinate system in which a 2 × 2 square matrix has nonzero off-diagonal elements. We do this by performing a similarity transformation on the matrix . The coordinate systems and are related by the similarity transformation matrix such that

The two coordinate systems are shown in **Fig. 1**.

An angle (*θ*) is associated with the transformation in **Eq. 47** by writing the 2D coordinate transformation as

The coordinate systems and are related by the counterclockwise rotation shown in **Fig. 1**. We have an aligned coordinate system with the principal axes of the permeability tensor. The diagonal tensor in the coordinate system has the form

....................(49)

where *k*_{max} is the maximum permeability in the direction *y*_{1} and *k*_{T} is the permeability that is transverse to k max in the direction *y*_{2}. We want to know how the elements of the permeability tensor change if we transform to the different coordinate system .

The relationship between the elements of and is

If we multiply **Eq. 50** from the left by ^{-1} and from the right by , we obtain

We find the elements of ^{-1} by solving

The result is

where is the transpose of . Substituting **Eqs. 51** and **54** into **52** gives

We can use **Eq. 55** to calculate the elements of for any rotation angle *θ*. If the permeability is isotropic, we have *k*_{max} = *k*_{T} = *k*_{iso} and **Eq. 55** simplifies to the form

In the special case of isotropic permeability, the orientation of the coordinate system does not affect the values of the elements of the permeability tensor.

**Fig. 2** shows the results for an anisotropic case in which *k*_{max} = 200 md and *k*_{T} = 50 md. The values of elements *k*_{11}, *k*_{12} and *k*_{22} of are presented in the figure. The off-diagonal terms satisfy the equality *k*_{12} = *k*_{21} for a symmetric matrix given in **Eq. 16**, so it is sufficient to show only *k*_{12}. The values of the diagonal elements change most as *θ* approaches 45°, and the values of the off-diagonal elements are greatest at *θ* = 45°. A rotation of 90° recovers a diagonal permeability tensor, but k max is now aligned along *x*_{2}, and *k*_{T} is aligned along *x*_{1}.

### Gridding a channel sand

The ideas discussed are now considered in the context of a realistic application. Our problem is to find a coordinate system that lets us accurately model fluid flow in a channel sand with the two permeability regions in **Fig. 3**.

We can highlight important features of the relationship between grid orientation and the assumption of diagonalized permeability by assuming the permeability in each region is anisotropic with a maximum permeability k max and a permeability *k*_{T} that is transverse to the direction of *k*_{max}. The diagonalized, anisotropic permeability tensor ′ in the *y*_{1} - *y*_{2} plane of a channel sand is the matrix

and Darcy’s law for flow in the *y*_{1} - *y*_{2} plane is

Consider two cases. In Case A, permeability is homogeneous and anisotropic, and in Case B, permeability is heterogeneous and anisotropic. The two cases are illustrated in **Figs. 4 and 5**.

**Fig. 6** shows two possible coordinate systems for orienting the grid. Coordinate system is more closely aligned with the spatial orientation of Region I than coordinate system , while coordinate system is more closely aligned with the spatial orientation of Region II than coordinate system . The coordinate system is obtained by rotating the coordinate system through an angle *θ* as in **Fig. 1**.

We consider four grid orientations for each case: (1) grid in Regions I and II; (2) grid in Region I and grid in Region II; (3) grid in Region I and grid in Region II; and (4) grid in Regions I and II. The grid orientation cases allow us to consider the effect of different coordinate systems on the permeability tensor in each region. We assume in our analysis that the reservoir simulator is a typical simulator with a formulation of fluid flow equations that uses Darcy’s law with a diagonal permeability tensor. We also assume the simulator allows different grid orientations in different regions of the model; otherwise, grid orientation cases 2 and 3 are not feasible. Our analysis does not include multidimensional numerical dispersion,^{[5]} which can also affect the accuracy of flow calculations. Results of the analysis are summarized in **Table 1**.

An "ok" in the "Permeability Tensor" column in **Table 1** indicates that the diagonal permeability tensor is aligned with the grid. An "X" indicates that the magnitudes of the diagonal terms in the permeability tensor must be corrected with **Eq. 55**. An "ok" in the "Formulation" column in **Table 1** indicates that the formulation of the fluid flow equations is correct. An "X" indicates that the formulation of the fluid flow equations is incorrect because the formulation does not include off-diagonal terms in the permeability tensor. Based on the results in **Table 1**, we observe that the grid orientation in Case A.1 provides the most faithful representation of the permeability tensor in Case A, and the grid orientation in Case B.2 provides the most faithful representation of the permeability tensor in Case B.

## Nomenclature

## References

- ↑ Gupta, A., Penuela, G., and Avila, R. 2001. An Integrated Approach to the Determination of Permeability Tensors for Naturally Fractured Reservoirs. J Can Pet Technol 40 (12): 43. PETSOC-01-12-02. http://dx.doi.org/10.2118/01-12-02
- ↑ Settari, A., Walters, D.A., and Behie, G.A. 2001. Use of Coupled Reservoir and Geomechanical Modelling for Integrated Reservoir Analysis and Management. J Can Pet Technol 40 (12): 55. PETSOC-01-12-04. http://dx.doi.org/10.2118/01-12-04
- ↑ Young, L.C. 1999. Rigorous Treatment of Distorted Grids in 3D. Presented at the SPE Reservoir Simulation Symposium, Houston, 14-17 February. SPE-51899-MS. http://dx.doi.org/10.2118/51899-MS
- ↑ Fanchi, J.R. 2006. Math Refresher for Scientists and Engineers, third edition. New York: Wiley Interscience.
- ↑ Fanchi, J.R. 1983. Multidimentional Numerical Dispersion. SPE J. 23 (1): 143–151. SPE-9018-PA. http://dx.doi.org/10.2118/9018-PA

## Noteworthy papers in OnePetro

Ramey, H.J., Jr. Interference Analysis for Anisotropic Formations - A Case History (includes associated paper 6406 ). http://dx.doi.org/10.2118/5319-PA.

## Noteworthy books

Carslaw, H., & Jaeger, J. (1986). Conduction of heat in solids (2nd ed.). Oxford Oxforshire, New York: Clarendon Press. 85026963 Worldcat

Earlougher, R. C. 1977. Advances in Well Test Analysis, SPE Monograph Series Vol. 5, Society of Petroleum Engineers, Richardson, TX, 264 pp. SPE Bookstore or Worldcat

Hahn, D., & Özisik, M. (2012). Heat Conduction. (3rd). Hoboken, New Jersey: John Wiley & Sons, Inc. Worldcat

Kreyszig, E. O. (2011). Advanced Engineering Mathematics. (10th). New York, New York: John Wiley & Sons Inc. Online Resource or Worldcat

## External links

Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro