Darcy’s law, The law governs one-dimensional flow of an incompressible homogeneous fluid in a homogeneous, isotropic porous medium. Our objective in this analysis  is to extend this law to other cases.

EXTENSIONS OF DARCY’S LAW

 The experimentally derived equation of motion, or Darcy’s law

 

 

(with J_s denoting the hydraulic gradient in the s-direction) is limited to:

  One-dimensional flow.

  Homogeneous isotropic porous medium.

  Constant density fluid.

  Nondeformable solid matrix.

  Flow at low Reynolds numbers.

 

 However, on the basis of the theoretical analysis e.g., Bear and Bachmat, 1990), that shows Darcy’s law to be a simplified form of the averaged momentum balance equation, the law presented above can be generalized to

 three-dimensional flow, 

 anisotropic

 inhomogeneous porous media, and 

 fluids of variable density.

 

 

 The objective of this analysis is to extend Darcy’s law (actually, the motion equation, as Darcy never referred to other cases) beyond the above constraints.

 

Three-dimensional flow

The specific discharge is a vector, q:

 

 

A vector is a quantity that has a magnitude and a direction in space. It can be

 

represented by its three components, say

 

q_x, q_y and q_z

 

in the x,y,z Cartesian coordinate system.

Note the difference between the symbols used for denoting the scalars x,y,z, and q, (i.e., regular italics) and those used for the vectors x, and q (i.e., bold-face italics).

 The gradient of a scalar function, say of h = h (x, y, z), denoted as “grad h“, or “” (to be read as “gradient (of) h“), is a vector that indicates at any given point, x, y, z in space, the magnitude and direction of the steepest ascent of the function. This is an extension of the concept of a derivative in one dimension. We recall that a derivative gives the slope of the function at the considered point. A positive derivative means that the function is rising at the considered point.

 

The hydraulic gradient is a vector equal to the negative of the gradient vector (i.e., the same magnitude, in the opposite direction). Thus, the hydraulic gradient at a point is a vector that indicates, in magnitude and direction, the steepest ascent, or the steepest rate of decline of the considered function at that point in the considered spatial domain.

 

 The following figure shows the contours of the function h = h(x,y) in a two-dimensional space—a plane. We see three contours: h₁, h₂, and h₃.

 

 

Fig. 1: Definition sketch for a gradient and a hydraulic gradient.

 

 

We see curves of equal values of the scalar h. The two-dimensional figure describes a sloping curved surface in space. It gives the value of h at any point x,y.

 

At any point in the x,y-plane, we may ask

 

“what is the steepest slope of the considered surface?”

 

 The direction of the slope at each point is given by the direction of the normal to the h-contours at that point. Along this normal, the gradient of h = h(x,y) is a vector pointing from low to high values of h.

 

 The hydraulic gradient at a point is a vector that has the same magnitude as the gradient, but points in the opposite direction.

 In a two-dimensional domain in the horizontal x,y-plane, we can show the gradient in a graphical form. Note the curves of equal piezometric head (like the elevation contours of a given hill):

 

 

 

 

Note that the steepest ascent means the shortest path from any contour to the higher one. The gradient is the ratio between the difference in piezometric head (1 m in the figure) and the shortest path between the contours, at the considered point.

 

 

 The following sketch presents the same idea in a different form. Note the vectors (in red) of the gradient of h, and their components in the x and y directions. The sketch shows a streamline that is everywhere normal to the contours. Note that the two observation wells are not along a streamline. .

 

 

 We recall that:

 

 The velocity (and the specific discharge) vector is everywhere tangential to the streamline, and

 In an isotropic porous medium, the streamlines are everywhere normal to the equipotential surface.

 

Here, an “equipotential surface” means “a surface at every point of

 

which the piezometric head is the same”.

 

 

 In the following sketch we note the gradient of the piezometric head at two points along a streamline:

 

 

The arrows indicate the specific discharge, which is tangent to the streamline, but they also indicate the steepest descent in the piezometric head at the considered points.

 

 

Again, on this sketch, we note that the specific discharge vector is tangential to the streamline, and that the local hydraulic gradient (As an average between the two points) is expressed by Dh/L.

 

 

 When flow takes place through a homogeneous isotropic porous medium domain, the coefficient K is a constant scalar, and the following three scalar equations maybe written as a consequence of Eqn. (1):

 

 

where x, y, and z indicate the axes of a Cartesian coordinate system, and

 

 

are the components of the hydraulic gradient vector.

 

 Altogether, the two sets of ordered numbers

 

q_x,q_yq_z, and , ,

 

are components of the vectors of specific discharge, q, and of the hydraulic gradient, J, respectively Note that the vectors qand J are independent of the selected coordinate system, while the magnitudes of their components do o the selected coordinates. Note that in the previous lecture w used the symbol J to denote the hydraulic gradient (denoted here by J).

 

 

 Hence, as an extension of the one-dimensional Darcy’s law, the three-dimensional Darcy’s law takes the form

 

 

or, in terms of the piezometric head, h,

 

 

 Note that ““, like “grad h“, is read as `grad (or gradient of) h‘. In mathematics, a function of space and time whose gradient gives the velocity of a continuum is called a potential. When the velocity is proportional to the gradient of some function, the latter is referred to as a quasi-potential. In view of Eqn. (3), the piezometric head, h, is a “quasi-potential”, although many groundwater hydrologists refer to it as `groundwater potential‘.

 

 

 The velocity vector, V, is given by

 

 

At a point, is a vector normal to the equipotential surface (or curve in two dimensions), i.e., a surface (or a curve) of constant h at that point. It points in the direction of steepest ascent in piezometric head. Hence, the vector = – indicates the direction and magnitude of steepest descent in the piezometric head. Note that we are not using any special symbol to indicate that is a vector, because by definition, the gradient of a scalar is a vector.

 

Try now to solve the following problems:

 

Problem 1:

 

Given the following observations of the piezometric heads in three observation wells (= piezometers):

 

Well	A	B	C
Coordinate x	0	300 m	0
Coordinate y	0	0	200 m
Piezometric head (m)	+ 10	+11.5	+8.4

 

Assume that the wells penetrate a homogeneous, isotropic, confined aquifer of constant thickness B = 20 m, porosity f = 0.2, and hydraulic conductivity K = 15 m/d, and that the piezometric surface between the wells can be approximated as a plane.

 

Required: determine (a) the hydraulic gradient (magnitude and direction), (b) the total discharge in the aquifer, per unit width, and (c) the velocity of the water at point P (100m, 100m).

 

ANSWER

 

Compressible fluid

 

 In a compressible, homogeneous (i.e., of constant concentration of dissolved solids) fluid, under isothermal conditions, the density is a function of the pressure, p, only, viz., r = r(p). For such a fluid, the motion equation (we still refer to it as Darcy’s law) is written in terms of Hubbert’s potential, h*, defined (on Darcy’s law) in the form

 

 

Thus, Darcy’s law takes the form

 

 

Inhomogeneous porous medium

 

As stated before, a domain is said to be homogeneous with respect to a porous medium coefficient if the latter is the same at all points in the porous medium domain. Otherwise, it is said to be heterogeneous or inhomogeneous. Unless otherwise specified, we shall use the term `homogeneous’ when the permeability of a considered porous medium domain is uniformly distributed, and so are other porous medium properties that are consequences of the structure of the solid matrix.

 

 Because these coefficients express behavior at a point within a domain, Darcy’s law in the form of Eqns. (1a) and (2) presented above remain valid also for three-dimensional flow through an isotropic heterogeneous porous medium domain, with the hydraulic conductivity at a considered point described by

 

K = K (x, y,z ).

 

Thus, Darcy’s law in a inhomogeneous domain is:

 

q = -K (x, y, z) .

 

 

Anisotropic porous medium

 

 

 If the permeability at a given point is independent of direction, the porous medium at that point is said to be isotropic. Otherwise, the porous medium is called anisotropic with respect to permeability

 For an anisotropic porous medium, Darcy’s law takes the form

 

 

 This equation expresses a linear relationship between the x, y, z components q_x, q_y, q_z of the specific discharge vector, q, and the components , J_x, J_y, J_zof the hydraulic gradient, vector J. These two sets of components are related to each other by nine coefficients, K_xx, K_xy,…,K_zz. Together they represent the hydraulic conductivity of the porous medium. In a heterogeneous medium, each of these nine coefficients may vary in space.

 

 The component K_ij, with i,j = 1,2,3, or x,y,z, may be interpreted as the contribution to the specific discharge in the ith direction, q_i, produced by a unit component of the hydraulic gradient in the jth direction, J_j. The total specific discharge in any direction is the sum of the partial specific discharges in that direction caused by J_x,, J_y,, andJ_z.

 

 The nine coefficients appearing in Eqn.(5) are components of the second rank tensor of hydraulic conductivity of an anisotropic porous medium,K..

 

What is a TENSOR????

 

 You have to know when a symbol denotes a vector, and when it denotes a tensor. In vector notation, we use roman boldface symbols (e.g., K for hydraulic conductivity), to denote second rank tensors, and boldface italics for vectors (e.g., q for specific discharge). The scalar components of both vector and tensor quantities are denoted by regular italics fonts, with subscripts that denote coordinates.

 

 A detailed discussion on the nature and use of second rank tensors is beyond the scope of this course. I refer you to texts on tensor analysis, such as Morse and Feshbach (1953) and Aris (1962).

 

 NEVERTHELESS, let me say a few words. You may have encountered the concept of a stress (especially if you are an engineer). If I ask you “what is a stress?” you may want to answer: “force per unit area.”. However, how can you divide a vector (i.e., force) by another vector (the area, which has a magnitude and direction)? Such quotient IS NOT DEFINED!

 

Instead, we may say that a FORCE VECTOR (F) is obtained by the product of stress by the AREA VECTOR (A). The stress (s) is then said to be a 2nd rank symmetric tensor:

F = A. s .

 

A vector (three components in a 3-d space) may be considered as a first rank tensor, while a scalar (one component) is a zero-rank tensor. The 2nd rank tensor, in a 3-d space has (or is described by) 9 components.

 

 

Another (more mathematical) way to explain what is a tensor, is to refer to its behavior under rotation of the coordinate (x,y,z) system. You are probably familiar with what happens to a vector (say, a force vector) at a point in space, as we rotate the coordinates. Obviously, the force (vector) itself is not affected by the fact that we rotate the coordinate system. Only the COMPONENTS of the vector vary as a consequence of this rotation. In mathematics, a certain rule expresses the relation between the magnitude of vector component before aafter rotation of the coordinate system. We may now turn this statement around and state that any three numbers that, under rotation of coordinates, behave according to this mathematical rule are components of a vector. Similarly, any 9 numbers that, under rotation of the coordinate system, behave according to the rule that defines the changes in magnitude of components of a second rank tensor, are components of a second rank tensor.

 

One way or another, the hydraulic conductivity in an anisotropic porous medium is a second rank tensor.

 

 The components of the tensor K in a three-dimensional space, can be written in the matrix form

 

 

and in a two-dimensional space, as

 

 

 

 The hydraulic conductivity tensor is symmetric, that is

 

 

This means that actually only six distinct components are needed to fully define the hydraulic conductivity in a three-dimensional domain (and only three in a two-dimensional one). Furthermore, the coefficients are non-negative.

 

 Equations (5) may be written in several compact forms. For example, they may be written in the compact vector form

 

 

i.e., expressing (the vector) q as a scalar (or dot) product of the (tensor of) hydraulic conductivity, K, and the hydraulic gradient (vector), J (= – ).

 

In indicial notation for an orthogonal coordinate system, Eqn. (8) can be rewritten as

 

 

where the subscripts i and j indicate x_i and x_j, respectively, with x₁ = x, x₂ = y, and x₃ = z.

 

 

In writing Eqn. (9), we have made use of Einstein’s summation convention (or double index convention) according to which a subscript (or superscript) repeated twice and only twice in any product or quotient of factors is summed over the entire range of values of that subscript (or superscript), i.e., i, j= 1, 2, 3 and 1, 2 for three- and two-dimensional spaces, respectively.

 

 Thus, in three dimensions,

 

 

 In this course, we shall often write equations in their compact vector form. Although only scalars are used when actually making calculations (and we recall that components of vectors are scalars), the advantage of writing the vector equations (in addition to the saving in the writing that is involved) is that given a problem, the equations can easily be rewritten in the coordinate system that is suitable for that problem.

 

 

 Although the hydraulic conductivity tensor, K, at a point within an anisotropic porous medium is independent of the coordinate system used, the magnitude of each K_ij-component does depend on the chosen coordinate system. Texts on tensor analysis give the rules for transforming these components from one coordinate system to another.

 

 Principal Directions

 

These texts also prove that it is always possible to find three mutually orthogonal directions in space such that when these directions are chosen as the coordinate system for expressing the components, we find that

 

K_ij=0 for all ij, and

 

K_ij 0 for i = j.

 

 These directions in space are called principal directions of the hydraulic conductivity tensor of the anisotropic porous medium. In a heterogeneous porous medium domain, the principal directions may vary from point to point.

 

 

 When the principal directions are aligned with the selected coordinate system, Eqns. (6) and (7) may be represented in matrix form as:

 

 

respectively, so that Eqn. (5) reduces to

 

q_x = K_xxJ_x, q_y = K_yyJ_y, q_y = K_zzJ_z(10)

 

 

 Because the hydraulic conductivity is related to the permeability by the scalar factor rg/m, the permeability of an anisotropic porous medium is also a second rank tensor. In fact, it is the tensorial nature of the permeability that determines the tensorial nature of the hydraulic conductivity.

 

 Let us solve some problems associated with anisotropy.

 

 

Problem:

 

Let the aquifer of an earlier problem be anisotropic:

 

Given the following observations of the piezometric heads in three observation wells (= piezometers):

 

Well	A	B	C
Coordinate x	0	300 m	0
Coordinate y	0	0	200 m
Piezometric head (m)	+ 10	+11.5	+8.4

 

Assume: The wells penetrate a homogeneous, anisotropic, confined aquifer of constant thickness B = 20 m, porosity f = 0.2, and hydraulic conductivity (in m/d)

 

 

 

and that the piezometric surface between the wells can be approximated as a plane.

 

Required: determine (a) the specific discharge, q, (b) the angle, a, between qand J.

 

Solve the problem ON YOUR OWN, and then compare with my ANSWER.

 

Generalized Darcy law

 

 We started the discussion on Darcy’s law or the equation of motion of a Newtonian fluid phase in a porous medium, from theoretical considerations. A review of this approach is beyond the scope of this book (see, for example, Bear, 1972; Bowen, 1976; Hassanizadeh and Gray, 1980; Bear and Bachmat, 1991). Although a number of different approaches have been employed, most of these investigations start by recognizing that the motion equation is a momentum balance equation of the fluid phase. They derive the (macroscopic) motion equation for a Newtonian fluid that moves in the void space of a porous medium, regarded as a continuum, by taking an average of the microscopic momentum balance equation of that fluid, known as the Navier-Stokes equation, over the volume of the phase within a representative elementary volume (REV) of the porous medium.

 

 By introducing a number of simplifying assumptions and, especially, by assuming that

 

 

• the inertial effects are negligible relative to the viscous ones, and
• the effect of gradients in specific discharge are negligible in comparison with the drag produced at the fluid-solid interface,

 

 

they obtained a simplified form of the averaged momentum balance equation.

 

 For the case of a single Newtonian fluid that occupies the entire void space (saturated flow), they obtain

 

 

where V, p, r and, m denote the (intrinsic phase) average velocity, pressure, density, and viscosity of the fluid, respectively, V_s denotes the (intrinsic phase) average velocity of the solid, z is the elevation, and k is the permeability tensor, which reduces to a scalar, k, in an isotropic porous medium. Note that at the macroscopic level, Vis what we usually refer to as the fluid’s velocity.

 

 

 In indicial notation and Cartesian coordinates, Eqn. (11) takes the form

 

 

Eqn (12) represents what happens in the ith direction.

 

QUESTION:

 

(a) How many terms do you see on the r.h.s?

 

(b) How many equations are represented by Eqn. (12) in a three dimensional space?

 

ANSWER.

 

Equation (12) is the general motion equa (we may still refer to it as Darcy’s law) for saturated flow of a single fluid in an anisotropic and heterogeneous porous medium, i.e., where k = k (x, y, z).

 

 

 Note the new feature: Darcy’s law gives the velocity of the fluid relative to that of the solid. Yes, the solid may be displaced, e.g., in the case of land subsidence, or consolidation.

 

 

 A comparison of Eqn. (11) with the various forms of Darcy’s law appearing earlier in this lecture, leads to the conclusion that the velocity in the latter is also a mass weighted one. This has no practical significance until we consider a multi-component fluid with spatial variations in component concentrations that produce density variations and fluxes by molecular diffusion. So, for the time being, you may disregard this comment. I’ll return to it in due time.

 

 

 In Eqn. (11), the term represents the macroscopic driving force per unit volume of fluid, due to a pressure gradient and to gravity.

 

 

The above figure shows the force balance at a point in a porous medium domain. You see the force resulting from the pressure gradient and that due to gravity. You see the resultant which is the vector sum of the two forces.

 

This resultant produces the flow as described by Darcy’s law.

 

This force is balanced by the drag, or resistance to flow at the solid-fluid interface. The latter is expressed by

 

 

 Written in this form, Darcy’s law states that the so-called Stokes drag is proportional to the fluid velocity relative to that of the soil skeleton, proportional to the fluid viscosity, and inversely proportional to the macroscopic coefficient k/f, which represents (at the macroscopic level) the effect of the microscopic configuration of the void space of the porous medium on the saturated flow.

 

 

Comment: the following two paragraphs are of interest, but if you have difficulties in understanding them, just skip them and move on.

 

 We wish to emphasize again that the motion equation for a fluid phase, e.g., Eqn. (11), is nothing but a simplified form of the macroscopic momentum balance equation for that phase. It expresses the exchange of momentum between the fluid and the solid phases.

 

It is also worth noting that although some forms of the motion equations (e.g., Eqn. (8)) look like a linear relationship that expresses a diffusive flux, Darcy’s law is definitely not an expression of this kind. The (phenomenological) laws of Fick, Fourier, and Newton express the diffusive fluxes of mass of a component, heat, and momentum, respectively. These laws exist at the microscopic level, and can be averaged to yield their macroscopic counterparts. Darcy’s law cannot be obtained by averaging a microscopic (phenomenological) diffusive flux law. It is obtained as an average of the microscopic momentum balance equation, and then simplified to by neglecting certain phenomena. The resulting law expresses the advective flux of the fluid’s mass, and not a diffusive flux.

 

 

 Introducing the definition of specific discharge relative to the solid

 

 

Eqn. (13) becomes

 

 

 Note in the above equation, as in earlier equations the “dot”. It is important. You should not omit it. It indicates a “dot (= scalar) product of two vectors”.

 

 In this equation, the fluid’s density, r, may depend on pressure, concentration of components, and temperature.

 

 For a compressible fluid, r = r(p), the right hand side of Eqn. (14) reduces to –K h*, where K=krg/m and h*, is Hubbert’s potential defined by Eqn. (5) of TOPIC B, LECTURE 1. For a constant r, the r.h.s. reduces to –K .

 

When V_s=0 (or approximately so, i.e., a stationary non-deformable solid), by replacing fV by q, we obtain

 

 

 When the x, y, and z coordinate axes are principal directions of the permeability of an anisotropic porous medium, and the direction of z is parallel to and in the sense opposite to that of gravity action, this single vector equation may be rewritten in a Cartesian coordinate system in the form of three scalar equations

 

 

 

Non-Darcian motion equations

 

This subject was discussed in detail in TOPIC B, LECTURE 1. Let us add a few more details.

 

Darcy-Forchheimer equation

 

 

 In the range of validity of Darcy’s law, i.e., Re < 1 (or even <10), the viscous forces that resist the flow are predominant. As the flow velocity increases, a gradual transition is observed from (macroscopically) laminar flow, where viscous forces are predominant, to still laminar flow, but with inertial forces gradually taking over. Often, the value of Re = 100 is mentioned as the upper limit of this transition region in which Darcy’s linear law is no more valid. The reason for this deviation from the linear law is that at the microscopic level, as velocities increase, local separations of the flow from the solid walls of the solid matrix occur in increasing number of places at which the flow curves or diverges. Local vortices and countercurrent flow regimes are caused by inertial forces at the microscopic level along portions of the solid.

 

 For flow in this range, the motion equation, known as Darcy–Forchheimer equation, is used. For a rigid porous medium, neglecting the effects of inertia at the macroscopic level, it takes the form

 

 

where q=|q|, n is the kinematic viscosity of the fluid, and the bwith components b_ij, is a tensor coefficient that is related to the configuration of the void space. At low Re, the second term on the r.h.s., which expresses the average of the microscopic inertial effects, becomes negligible. Then the equation reduces to the linear Darcy law.

 

 

Brinkman’s equation

 

 

When all inertial effects are negligible, but we do not neglect the dissipation of energy within the fluid, as fluid layers move at different velocities, the motion equation takes the form known as Brinkman’s equation (Brinkman, 1948)

 

 

where T*_pq denotes the pq-component of the tortuosity (2nd rank) tensor, and all other symbols have already been explained earlier.

 

 

Rapid velocity variations

 

 

 When, especially at the onset of flow and in oscillatory flows, local acceleration may not be neglected, but we do neglect the advective acceleration and the internal friction within the fluid, the motion equation takes the form

 

 

 

We note the effect of acceleration, .

 

 Although we have introduced these two cases here for the sake of completeness, conditions that justify their application are seldom, if ever, encountered in problems of flow and contaminant transport in the subsurface considered in this course. We do not intend to make use of these flux equations in our course.

 

Let’s summarize:

Altogether, given a problem of flow in a saturated porous medium domain, the appropriate form of the motion equation depends on the answers to the following questions, which are posed in the development of the conceptual model:

 

 Are inertial forces negligible (i.e., Re < 1 to 10), so that (the linear) Darcy’s law can be applied?

 Is the porous medium homogeneous, or heterogeneous, isotropic, or anisotropic?

 Is the fluid’s density constant, or is r = r(p), orr = r(p, c, T)?

 

 

Once we have answered these questions, there should be no difficulty in selecting the appropriate form of the motion equation from among those presented above.

 

 

 

 

 

 

 

 

 

 

 

ANSWER

Answer to Problem 1(a)

 

A. The problem can be solved in three ways.

 

 (1) We start by drawing the equipotentials on the map, using linear interpolation. We always use linear interpolation unless we have a specific reason to assume otherwise. What we actually say is that the three points define a planar piezometric surface. The hydraulic gradient (magnitude and direction) is then the slope of this plane. It is determined by measurements on the map. Because the three values of the piezometric head constitute a plane with a constant slope.

 

 

The above figure shows the interpolated contours. The thick arrow, perpendicular (why?) to the contours indicates the direction of the vector of hydraulic gradient, J (note that it is from high to low piezometric heads). Its magnitude can easily be determined from this map, as the ratio Dh/L (i.e., ratio between the Dh between any two contours, to the distance between the contours.

 

Note that we have selected an arbitrary point within the considered ABC domain to indicate this arrow. Does the selection of the point (say where the arrow starts) affect the magnitude of the gradient?)

 

 (2) A second option is for the more mathematically minded.

 

The shape of the piezometric surface is a plane (what else can we choose as a surface passing through three given points, in the absence of any other information?)

 

h = h(x, y) = ax + by + c,

 

in which a,b,c are constants. By substituting the given data into this equation for each of the three points, we obtain:

 

10 = a x 0+b x 0+c, hence, c=10.

 

11.5 = a x 300+b x 0+10, hence, a=5 x 10^-3.

 

8.4=a x 0+b x 200+10, hence, b=-8 x 10^-3.

 

The equation of the piezometric surface is, therefore,

 

h =5 x10^-3 x-8 x 10^-3y+10.

 

We can now determine the hydraulic gradient and its components from the definition:

 

 

The third approach is similar to the first, but based on the given data: Because the piezometric surface is a plane, with the same slope everywhere, we can determine the slope directly from the figure, or from the given data:

 

 

(b) The total discharge:

Q’ = K B W = (15m/d) x (20m) x (1 m) x (9.43 x10^-3) = 2.8m3/m/d.

(c) The location of the point does not play any role, as the piezometric surface is a plane (i.e., the same slope everywhere), and the aquifer is homogeneous. Hence:

 

ANSWER

 

In an earlier problem, we have already found that in the case considered here:

 

Hence:

 

 

The specific discharge vector is:

where 1x and 1y are unit vectors in the x- and y-directions, respectively.

The absolute value is:      q = 72.14m/d

 

(b)

In the same earlier problem, we obtained , as the angle that the vector J makes with the +x axis. The angle that the flux vector, calculated in (a), makes with the same axis is:

Hence:

 

ANSWER

(a) On the right hand side, for the same i, we have a sum of three terms: for j = 1 plus a term for j = 2, and a term for j = 3.

(b) In Eqn (12) we keep the value of i.The value of j varies: j = 1j = 2, and j = 3.

However, i may take on the values 1, 2, 3. Hence the equation really represents three equations.

REFERENCE BOOKS:

*Modeling Groundwater Flow and Contaminant Transport By Jacob Bear, A. H.-D. Cheng (PG.146)

* Fractals in Petroleum Geology and Earth Processes.,edited by Christopher Cramer Barton, P. R. La Pointe (PG.188)

http://www.me.ust.hk/~mezhao/pdf/19.PDF

*Advanced Petroleum Reservoir Simulation, By Rafiq Islam, S.H. Moussavizadegan, Shabbir Mustafiz, Jamal H. Abou-Kassem

http://echo2.epfl.ch/VICAIRE/mod_3/chapt_5/main.htm

http://www.petrobjects.com/downloads/Petroleum%20Simulation%20Papers%20&%20Articles/Design%20&%20Implementation%20of%20a%20Numerical%20Model%20for%20Simulating%20Petroleum%20Reservoirs.pdf

http://www.slb.com/~/media/Files/resources/oilfield_review/ors96/sum96/06961627.pdf

 

 

 

 

 

 

 

 

 

What is the difference between numerical solution and approximate solution. Also when we can say that the method is numerical or approximate method?

Comments

All numerical solutions are approximate.

All approximate solutions are numerical.

When trying to make a model of the actual physical reservoir, all models are approximations, and all models have to be solved numerically.

A numerical solution is when you take the physical equations and solve them using numerical methods or algorithms, i.e use approximation for differential operators assuming small increments. dx/dt = (x1-x2)/dt

An approximate solution is when you solve the physical equations exactly but use approximate values to get the final answer i.e solve the differential equation and then use approximate values. i.e solve dx/dt=2 therefore x=2t +C and use approximate values for x and/or t.

numerical solution are usually applied to highly non-linear equation, e.g complex flow partial differential equation. It might involve the application of finite difference approximation using Taylors series or other discretization algorithm.

Approximate solution applies when you try to obtain a direct solution to a complex polynomial or equation which may involve unequal fractional powers. I always use approximate solution anytime I am dealing with Juhasz or Waxman-Smith water saturation equations.

There can be numerical and analytical solutions. Analitical solutions can be both precise and approximate. Numerical solutions are always approximate.

Problem of analitical solutions is that they can only be applied to very simple problems (i.e. simple geometry of the well and reservoir). So even if the analitical solution is precise the answer you get at the end of the day is approximate because in most cases before using analytical solytion you had to simplify the problem.

When it comes to the numerical solutions in flow simulation there are two subtypes of the models:

“Full physics” models – they pretend to be as precise as it can be. I.e. in flow simulators we try to account for all physical effects we know how to account for. But even if this case it is hard to say that any particular model does “everything possible”. For example you can always refine the grid i a hope to improve solution.As well Full Physics models are known for long simulation times – that may be a problem for uncertainty modeling and assisted history matching in particular.

On the other end there are reduced physics models when we neglect some effects (i.e. gravity, heterogeniety PVT etc.). One simple way of getting “reduced physics model” is coarsening “reduced physics model”. Those models can be very fast but if chosen incorrectly may introduce unrealistic results (especially when soing uncertainty quantification and assisted history matching).

Numerical and approximate solutions are based on the same physics of developing the mathematical equations needed to solve a problem. A coarse model of a reservoir can be considered as an approximate solution (approximate analysis) where a fine gridded model for the same reservoir is considered as a numerical solution (numerical analysis). However, mathematical equations used in both cases are still approximate as they are developed based on assumptions that not necessarily representing the real case e.g. modelling of homogeneous and heterogeneous reservoirs are based on the same mathematical models and amount and quality of structure and fluid data i.e. both are approximate. Reservoirs are systems that are mathematically mimicked based on the available data that by all means will always result in approximate solutions regardless of the amount and accuracy of data used to model those reservoirs.

I think we have may be in danger of making it too complicated.

We have a physical reservoir. We have a mathematical model of the reservoir.

The mathematical model, normally encoded in software, is an approximation to the physical reservoir.

We use the mathematical model to perform studies, and hope the results of these studies are close enough to the physical model to make sound decisions.

We can solve the mathematical model in many ways. If it is a very simple model, we can solve it analytically. In most cases, we have to solve it numerically.

By analytically, I mean you can write out all the equations and the solution can be expressed as a mathematical equation which can (in principle) be solved exactly. Plug the numbers into your equation and you have the solution.

But even here, you may have an analytical solution, but it may be expressed with complex functions, or even simple functions like sin, cos, exp, Bessel etc. and these are still solved using approximate iterative methods, so it is still numerical even if it appears to be analytical.

The analytical solution is an exact solution to the model, but the model is an approximation to the physical reservoir.

The numerical solution is an approximate solution to the model, and the model is also an approximation to the physical reservoir.

In most cases the important difference is not between analytical or numerical solutions (provided care is taken in the numerical solution, which it is in all commercial reservoir simulations), but between the mathematical model and the physical reservoir.

As an example, for a given SPE test problem, different commercial simulators (should) give very similar results, which shows that the numerical solution methods have high accuracy.

But nobody would say those SPE test problems bear much relationship to any known physical reservoir.

As an aside, I don’t know whether they do or not, but it would be comforting if there were models for which the analytical solution was known, and the commercial simulators modelled it numerically to reproduce the analytical solutions. The fact that they give similar results to the same model is not sufficient to demonstrate they give the correct results.

Uncertainty quantification is quite important to verify the numerical solutions by testing against known analytical problems (e.g. variance of simple correlated Gaussians in high dimensions), but I am not aware of any software which does this, hence the big differences we see in uncertainty quantification between commercial vendors.

If your numerical solution cannot replicate the analytical solution on simple problems, or you have not attempted to replicate, go away and come back when it does.

Is the question about the proxy models derived by exact numerical models. ?
If so the the answer could be different.

In order to make it simply, all models of any mathematical kind or physical models are not exacts to a reservoir as it is well known, but numerical or analytical equations of a mathematical model, if they are simple can give exact solutions to the models but not to the reservoirs(even if with feeding by stochastic reservoir data).

Using analytical, numerical or finite difference methods of solutions of complex system or equations are subject to the type of that system and equations(complexity of fluid and flow) of our mathematical model towards exact or approximate solution

What is the difference in petrel and eclipse simulation of the same reservoir model without upscaling it in Eclipse? (Considering Petrel as a simulator not as static modeling software)

Comments
Consider Petrel as a pre and post-processor of data from Eclipse (input data and results).

What does Petrel do when you invoke simulation module? Well it creates all files neede by eclipse, then submits these files to be run in Eclipse…and uploads the simulation results inro Petrel.
Let me be redundant…all data files are still run by Eclipse

are you talking about Petrel RE?… Static models are finer (in grids and details) than dynamic models (Eclipse), so there is a need of a scaling process, so you’ll have coarser grids and averaged properties that might enable you to run your model faster. Is that what you meant with your question?

No, its not about upscaling of cellular model but it was about the use of Petrel as a simulator and comparison with Eclipse.

Petrel is not a numerical simulator.

Petrel RE can create the files needed for ECLIPSE simulator and load the results back for analysis. When you run a numerical simulation from Petrel, you are invoking ECLIPSE to do the simulation.

It is the same thing other than Petrel RE is much easier to use. Maybe you lack some keywards, but you may find them on on Define Simulation, Advanced, Editor and than it would make your life much easier

Petrel is not a simulator, the simulator is Eclipse.. you can consider petrel as an integrated software where all the pre-processor and post-processor data is created and displayed.. Without Petrel you have to use other separated modules to create your data like FloGrid.. and modules to display your results like Floviz and Eclipse office

Eclipse is a Simulator and Petrel just makes static Geological Models and also Petrel convert the fine Model into Coarse by Upscaling method…an other difference is the grid used in Petrel is corner point which provide good flexibility for the shape of Geological Model while in Eclipse we can’t Simulate with this geometry of corner point cells as we are dealing with the flow rate which are controlled by Pressure gradient and by Darcy’s law so the gridpoint should be in the centre(midway btw the nodes)…..then we can get good flow behavior for simulation…..

If I understand Petrel is much used for the geological model and Petrel RE module who is under development by the owner Schlumberger provides a predefined idea of the dynamic model, but we can not do without use eclipse module is the most comprehensive dynamic simulation while knowing that the files are generated from petrel.

Maybe PETREL RE will replace Eclipse thereafter, it will have one complete software for static geological model and development model dynamics.

Petrel will never replace Eclipse as simulator … Schlumberger comme up with new simulator called Tempest which allows lot of flexibility in terms of grid shape …. Petrel is a pre and post precessor and Eclipse is the simulator …. The only good thing I found by running models through Petrel is the algorithm that generates transmissbities in the model … it make the model running bit faster and stable ……

Petrel will never replace ECLIPSE, Petrel is a post & pre processor for Eclipse simulation model.

The advantage now with Petrel is that you can run assisted history matching using the objective function & uncertainty & optimization….

Petre RE is not a simuator, I tis like an early easy interface to generate the eclipse files before simulation occurs ( in eclipse software), the developement strategy for example is a plus.

I just want to know if the upscaling in the Pertel RE is exactely the same as that one done by petrel ( if someone expereinced this before)

Petrel RE only as interface to bring all the input to the dek files and after that its still need enginge which is simulator called eclipse 100/300

We use Petrel RE for simulation due to the option of pre-post processing. Petrel RE has the same deck with the Eclipse Office and the simulation is run also with Eclipse. But the uncertainty and the optimization tools are very helpful if you think about the uncertainty analysis. If you are planning to have some intelligent well completions you can easily handle with well segmentation and installation of ICDs/ICVs in Petrel RE with the completion modules. Petrel RE has also the plug-ins like Mepo and Olyx. Therefore you can carry out a whole study from static modeling to the sensitivity analysis in one project. This will give you the option to have a better view of your results. It was also mentioned above almost every keyword is available also in Petrel RE but if you cannot find the keyword that you need you can manually add in the deck file (DATA file) or you can add/edit manually in Advanced sections. The only problem in Petrel RE is that if you are planning to run a Radial generic model, unfortunately you cannot easily create this geometry. You have to manipulate it

Tempest is nothing to do with Schlumberger, it is a simulator from Roxar. Maybe you are thinking of Intersect?

all above comments are valid about Petrel-RE, but I need to clarify the following:

1- Petrel-RE is never going to replace ECLIPSE for a simple reason that ECLIPSE is our simulation engine that solves the flow equations which Petrel is not going to do.
2- Our new Simulator is INTERSECT but is NOT a replacement for ECLIPSE. It’s used for very specific cases when having a highly fractured reservoir with millions of cells.
3- Schlumberger have nothing to do with Tempest 🙂

Tempest is a simlator of ROXAR Ltd. The logic and structure of data files of eclipse is very similar. However TEMPEST has less options for modeling than Eclipse but simulators results are converged

What do we really mean when we say that Petrel is a pre and post processor for eclipse simulation model?

Pre-processor because Petrel RE will help to create all files that Eclipse needs before submiting a run.

Post-processor because you may observe all results reported by Eclipse in either 3D grid or graphs, after Eclipse ends a run.

Can I run the files generated by a newer version of petrel (Petrel 2010) on an older version of Eclipse (Eclipse 2009) in the pre-processing stage. And if that’s possible can we do it the other way round, at the post-processing by exporting the older Eclipse simulation results into the newer Petrel??

The answer is yes, you can run DATA files generated with Petrel 2010 on an older version of Eclipse, provided that the keywords generated by Petrel are supported in your Eclipse release. The inverse is also affirmative. Eclipse results of any version can be loaded into Petrel.

i want to work on” Field development of gas condesate reservoir using compositional simulator.”as my project.i want yours help to work on simulator.which software in simulation should i use?

Petrel RE helps to you create your simulation model (that’s what we call pre processing), and then you can use Eclipse 300 (compositionnal) to run your model after that in Petrel RE you can analyse your results (line plots, 3D, 2D) and this is post processing.

I suggest to use PVTi or PVTsIM in order to calibrate your equation of state.

how can i load the Petrel model in Eclipse?

You cannot load the Petrel model in Eclipse unless you run the simulation. You can run the simulation in Petrel itself or you can export the Petrel keywords for your data deck. Those keywords are: GRID, PROP (porosity, permeability, water saturation), REGIONS ( FIPNUM, EQLNUM), SOLUTION (OWC, reference depth) etc. Once you run the simulation than you can see you 3-D model parameters in Eclipse.

Petrel RE is a visual aid for Eclipse code. You can run the simulation (waterflooding etc.) and visualize it. You can see the input data as plots (kr, pc etc) and make cross-plots (water cut, recovery etc.). Verry valuable tool.

you can also do the process in reverse. Create the model in Petrel and then export it to Eclipse, then re-import it for visualization.

I am working on ECLIPSE and Petrel.

ECLIPSE is a simulator.
Petrel is a pre and post processor tool. If you simulate a model in Petrel, it still uses ECLIPSE simulation engineer (or INTERSECT, new fast simulator). Therefore Petrel is NOT a simulator.
Petrel is easy for visualisation of the grid, simulation results etc.

In Russia that’s may be differences at less 5 % (after upscaling before hydrodinamic simulation) – thats percents are allowing, when Russian Central comissy for Production for oil fields is studing working simulations models of russian oil fields every 3 – 5 years

The question is without ‘UPSCALING’. That’s really a good question.

When you build the geostat model usually in Petrel, you try to honor the seismic, geological and petrophysical properties. To do that you need to build a high resolution Grid, which means with a large number of cells. Now let’s come to ‘upscaling’. Upscaleing means to reduce number of cells in order to run it in Eclipse. But that depends on the power of computer that you have. If you have a computer with 4-8 Gb RAM means your computer is not powerfull, so for a model above 100,000 cells you need to upscale. So my answer to your question is: it’s a matter to honor the model properties in compliance with seismic, geology and petrophysical. Upscaling usually means bad work.

Welcome to WordPress. This is your first post. Edit or delete it, then start blogging!