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 Js 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
qx, qy and qz
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: h1, h2, and h3.
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
qx,qyqz, 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).
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) .
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 qx, qy, qz of the specific discharge vector, q, and the components , Jx, Jy, Jzof the hydraulic gradient, vector J. These two sets of components are related to each other by nine coefficients, Kxx, Kxy,…,Kzz. 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 Kij, with i,j = 1,2,3, or x,y,z, may be interpreted as the contribution to the specific discharge in the ith direction, qi, produced by a unit component of the hydraulic gradient in the jth direction, Jj. The total specific discharge in any direction is the sum of the partial specific discharges in that direction caused by Jx,, Jy,, andJz.
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 xi and xj, respectively, with x1 = x, x2 = y, and x3 = 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 Kij-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
Kij=0 for all i
j, and
Kij 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
qx = KxxJx, qy = KyyJy, qy = KzzJz(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, Vs 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?
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 Vs=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 bij, 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.
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:
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:
(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